arXiv CSopen access

Average Gradient Outer Product in kernel regression provably recovers the central subspace for multi-index models

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

machine learning, deep learning, neural networks

Average Gradient Outer Product in kernel regression provably recovers the central subspace for multi-index models Libin Zhu∗

Damek Davis†

Dmitriy Drusvyatskiy‡

Maryam Fazel§

arXiv:2605.15082v1 [stat.ML] 14 May 2026

May 15, 2026

Abstract We study a prototypical situation when a learned predictor can discover useful low-dimensional structure in data, while using fewer samples than are needed for accurate prediction. Specifically, we consider the problem of recovering a multi-index polynomial f ∗ (x) = h(U x), with U ∈ Rr×d and r ≪ d, from finitely many data/label pairs. Importantly, the target function depends on input x only through the projection onto an unknown r-dimensional central subspace. The algorithm we analyze is appealingly simple: fit kernel ridge regression (KRR) to the data and compute the Average Gradient Outer Product (AGOP) from the fitted predictor. Our main results show that under reasonable assumptions the top r-dimensional eigenspace of AGOP provably recovers the central subspace, even in regimes when the prediction error remains large. Specifically, if the target function f ∗ has degree p∗ , ∗ it is known that n ≍ dp samples are necessary for KRR to achieve accurate prediction. In contrast, we show that if a low degree p component of f ∗ already carries all relevant directions for prediction, subspace recovery occurs in the much lower sample regime n ≍ dp+δ for any δ ∈ (0, 1). Our results thus demonstrate a separation between prediction and representation, and provide an explanation for why iterative kernel methods such as Recursive Feature Machines (RFM) can be sample-efficient in practice.

Introduction

Modern machine learning systems often succeed not only by fitting input-output relationships, but also by learning useful intermediate representations. This theme is central to contemporary representation learning: models are trained on high-dimensional data, yet the task-relevant structure is often believed to depend on a much smaller number of latent factors. A basic theoretical question is whether a trained predictor can reveal such low-dimensional structure with far less data than it needs for accurate prediction. In this paper, we study this question in a simple and tractable setting. We consider multi-index regression, in which the response y depends on input data x ∈ Rd only through a small number of linear measurements, y = f ∗ (x) + noise, f ∗ (x) = h(U x), (1) where U ∈ Rr×d has orthonormal rows with r ≪ d and h : Rr → R is an unknown link function. The row space of U , denoted row(U ), is called the central subspace. Recovering this subspace is valuable in its own right: it provides an interpretable representation of the data, enables dimension reduction, and reduces downstream learning to an r-dimensional problem. ∗ Department of Mathematics, University of Washington, Seattle, WA 98195; https://libinzhu.github.io/. Research of Zhu was supported by NSF TRIPODS II DMS-2023166. † Wharton Department of Statistics and Data Science, University of Pennsylvania, Philadelphia, PA 19104, USA; www.damekdavis.com. Research of Davis supported by an Alfred P. Sloan research fellowship and NSF DMS award 2047637 ‡ Department of Mathematics, U. Washington, Seattle, WA 98195; https://sites.google.com/view/dmitriy-drusvyatskiy. Research of Drusvyatskiy was supported by NSF DMS-2306322, NSF DMS-2023166, and AFOSR FA9550-24-1-0092 awards. § Department of Electrical & Computer Engineering, University of Washington, Seattle, WA 98195, and Amazon, Inc. https://people.ece.uw.edu/fazel_maryam/. Research supported by NSF TRIPODS II DMS-2023166, CCF 2007036, CCF 2212261, CCF 2312775 and the Moorthy Family professorship at UW.

Gradient information provides a direct population-level characterization of the central subspace. If f ∗ (x) = h(U x) is differentiable, then the population Average Gradient Outer Product (AGOP) takes the form E[∇f ∗ (x)∇f ∗ (x)⊤ ] = U ⊤ E[∇h(U x)∇h(U x)⊤ ]U, (2) and its range coincides with row(U ) in nondegenerate settings. This observation underlies averagederivative and outer-product-of-gradients methods in efficient dimension reduction. At the same time, this is only a population statement: the regression function is unknown, and naively estimating its gradient field in high dimension is typically no easier than solving the original nonparametric regression problem. Recent work on Recursive Feature Machines (RFM) uses AGOP-based updates for adaptive kernel machines, and has shown strong empirical performance on multi-index models and related learning tasks [8, 46, 47]. One aim of the present paper is to explain the feature-learning mechanism underlying this empirical behavior. To this end, we study the AGOP formed from a fitted kernel ridge regression (KRR) predictor, and investigate whether the gradients of the fitted predictor can reveal the central subspace even when the predictor has not yet learned the full target function. We show that the empirical AGOP constructed from a fitted KRR predictor provably recovers the central subspace in a high-dimensional regime. We focus on Boolean data drawn uniformly from {−1, 1}d , fit KRR with an inner-product kernel, and form the empirical AGOP matrix n

X c := 1 M ∇fˆ(x(i) ) (∇fˆ(x(i) ))⊤ ∈ Rd×d , n i=1

(3)

where fˆ is the KRR predictor and the gradients are ambient gradients evaluated at the training inputs x(i) . c concentrates around a low-degree population AGOP Our main results show that the empirical AGOP M target (Theorem 2), and under additional nondegeneracy and a weak subspace coherence condition [11], c consistently recovers the central subspace row(U ) (Corollary 4). the top-r eigenspace of M Subspace recovery before accurate prediction. The central message of this work is that recovering the central subspace from a fitted predictor can be statistically easier than achieving low test error. Specifically, in the sample regime n ≍ dp+δ with δ ∈ (0, 1), KRR primarily captures the degree-≤ p component of the target function. The prediction error can therefore remain large when a substantial portion of the signal lies in higher degree components. Our contribution is to show that this does not preclude accurate recovery of the central subspace: if the low-degree component already captures all relevant directions in the central subspace, then the gradients of the fitted predictor remain informative enough for the empirical AGOP to recover the subspace. Thus, the sample complexity of representation learning is governed by the first degree at which the low-degree gradient covariance contains all relevant directions, rather than by the degree needed to learn the full target function. This viewpoint suggests a natural two-stage approach. First, fit KRR and recover the central subspace row(U ) via the empirical AGOP. Second, project the input data onto the estimated subspace and learn the full link function in dimension r. Since r is small, the second stage is expected to learn higher-degree structure without incurring a d-dimensional sample-complexity cost. This viewpoint is also related to iterative kernel methods such as RFM [47] and IRKM [56], which aim to refine feature discovery over successive iterations. See the discussion in Section 4. Our experiments show that RFM identifies the central subspace well before its prediction error becomes small, consistent with our subspace-recovery results; subsequent iterations then improve prediction empirically. See Section 5 for numerical experiments. Contributions. Our contributions are threefold. First, we establish a subspace recovery guarantee for the empirical AGOP of the fitted KRR predictor in a polynomial high-dimensional regime. Second, we show a prediction–representation separation: the central subspace can be recovered accurately even in low-sample regimes where the prediction test error remains large. Third, we develop analytical tools for random matrix analysis with hypercube data, which may be of independent interest.

1.1

Notation

Indices and linear-algebra notation. We denote by N the set of natural numbers and write N+ := N \ {0}. For n ∈ N+ , let [n] := {1, . . . , n}. A vector λ ∈ Nd is called a multi-index, and we write |λ| := λ1 + · · · + λd . For x ∈ Rd and λ ∈ Nd , we use the standard monomial notation xλ := xλ1 1 · · · xλd d , Qd and we write Heλ (x) := i=1 Heλi (xi ), where Hek is the univariate probabilist’s Hermite polynomial of degree k. Pd Pd For u, v ∈ Rd , write ⟨u, v⟩ := i=1 ui vi for the Euclidean inner product and ∥v∥p := ( i=1 |vi |p )1/p for the ℓp -norm of v. For a matrix M , we denote its (i, j)-entry by Mi,j , its ith row by Mi,: , and its jth column by M:,j . We write diag(M ) for the vector of diagonal entries of M and Diag(v) for the diagonal matrix with diagonal v. For a matrix U ∈ Rr×d , the symbol row(U ) ⊆ Rd denotes its row space, and U⊥ ∈ R(d−rank(U ))×d denotes the orthonormal complement of row(U ). The symbols ∥M ∥op and ∥M ∥F denote the spectral and Frobenius norms of M , respectively. We write τd for the uniform probability measure on Hd := {−1, 1}d , and for a function f : Hd → R we R 2 2 define the usual square L2 -norm ∥f ∥L2 (τd ) := Hd f dτd . Differential notation. For a differentiable function f : Rd → R, we write ∇f (x) := (∂1 f (x), . . . , ∂d f (x)) for its Euclidean gradient. For any point on the hypercube x ∈ Hd , we regard x as a point of Rd and use the same notation ∇f (x) for the ambient gradient evaluated at x. Asymptotic notation. For nonnegative functions f and g, we write f ≲ g if there exists a numerical constant C > 0 such that f (x) ≤ Cg(x) for all x; similarly, we write f ≳ g to mean g ≲ f , and f ≍ g means both f ≲ g and f ≳ g. For functions f, g : Rd × N → R, the notation f = Od (g) means that there exist constants C > 0 and d0 ∈ N such that |f (x, d)| ≤ C|g(x, d)| for all x and all d ≥ d0 , while f = od (g) means that for every ϵ > 0 there exists dϵ ∈ N such that |f (x, d)| ≤ ϵ|g(x, d)| for all x and all d ≥ dϵ . In both cases, the comparison is uniform in x. Finally, for sequences of random variables fd and gd , we write fd = Od,P (gd ) if for every ϵ > 0 there exists Cϵ > 0 such that lim supd→∞ P(|fd | > Cϵ gd ) ≤ ϵ, and we write fd = od,P (gd ) if for every ϵ > 0, we have limd→∞ P(|fd | > ϵgd ) = 0.

Related work

Sufficient dimension reduction and multi-index models. Sufficient dimension reduction (SDR) studies regression settings in which the conditional distribution, or at least the conditional mean, of the response depends on a covariate x ∈ Rd only through a low-dimensional linear projection. Classical inverseregression methods include sliced inverse regression (SIR) and principal Hessian directions (pHd) [33, 34]; see also [13, 36].A complementary gradient-based line estimates dimension-reduction directions using average-derivative estimators [24, 45], outer products of gradients and MAVE-type local-linear procedures [52, 53], and structure-adaptive estimators [25]. Closely related are kernel SDR methods: kernel dimension reduction (KDR) characterizes sufficient directions through RKHS conditional covariance operators [20], whereas gradient-based KDR estimates directions from gradients of a kernel-estimated regression function [21]. Single- and multi-index models have also been studied extensively in semiparametric statistics and, more recently, under high-dimensional structural assumptions such as sparsity or smoothness [27, 35, 54]. Our work differs from this literature in both estimator and regime: we analyze the empirical AGOP of a trained KRR predictor and prove post-hoc subspace recovery in a polynomial high-dimensional regime. High-dimensional KRR and spectral structure. High-dimensional kernel methods have been analyzed through random-matrix and spectral descriptions of empirical kernel matrices and KRR risk [19, 44]. For inner-product kernels in polynomial high-dimensional regimes, analyses under spherical and hypercube covariate models show that KRR prediction is governed by low-degree polynomial components, leading to a degree-truncation phenomenon [22, 37]. We use this low-degree structure but study a different object: rather than characterizing the test error of the KRR predictor, we analyze the concentration and eigenspace structure of the empirical AGOP formed from its fitted gradients. 3

Feature learning in neural networks. Single- and multi-index models are standard testbeds for understanding feature learning in neural networks. For single-index targets, learnability under online SGD, gradient flow, and related shallow-network dynamics is often governed by Hermite coefficients, information exponents, or generative exponents [6, 10, 15, 32, 40, 51]. Related high-dimensional analyses show how first-layer adaptation can improve over fixed random features or kernel methods [7, 38]. For multi-index and sparse low-dimensional targets, recent work studies representation recovery [9, 16, 41, 42, 49]; staircase and leap-complexity phenomena [1, 2]; finite-step effects [17]; batch-size and batch-reuse effects [4, 5, 18]; and proportional-limit descriptions of feature learning [39]. The generative-leap framework of [14] establishes sharp thresholds for efficient hidden-subspace recovery in Gaussian multi-index models. In contrast, we study KRR with a fixed kernel—a setting closely related to the kernel/NTK regime of neural networks [28]—and show that the fitted gradients can contain sufficient information for post-hoc representation recovery, without weight evolution. Kernel-based feature learning and AGOP/EGOP methods. A growing literature shows that feature learning can also be realized through adaptive kernel methods. Recursive feature machines (RFM) [47] and iteratively reweighted kernel machines (IRKM) [56] alternate kernel regression with gradient- or AGOP-based metric updates, while related approaches optimize a linear map inside the kernel [26], use coordinate-wise kernel reweighting [48], or adapt local EGOP metrics, where EGOP is the population analogue of AGOP [30]. These works adapt the kernel or metric during training; in particular, Huang et al. [26] prove finite-sample excess-risk guarantees for HKRR, a representation-learning variant of KRR that optimizes a linear map inside the kernel, while we study the AGOP of a fitted ordinary fixed-kernel KRR predictor in the polynomial high-dimensional regime. Our closest comparison is [56]: both works use gradients of fitted kernel predictors under hypercube covariates, but [56] studies sparse and hierarchical targets via iterative metric learning, whereas we study general multi-index models and the eigenspace structure of the resulting AGOP.

Main Results

Throughout this section, we let the hypercube Hd = {−1, 1}d be equipped with the uniform measure τd , and all L2 -norms on Hd are taken with respect to τd . The main result has two steps. First, we show that the empirical AGOP of the fitted KRR predictor concentrates around a low-degree population target M≤p . Second, we show that if this low-degree target preserves the central subspace, then under a weak coherence condition, the top-r eigenspace of the empirical AGOP consistently recovers the central subspace row(U ).

3.1

Model and low-degree population targets

We consider the multi-index model f ∗ (x) = h(U x),

x ∈ Rd ,

(4)

where U ∈ Rr×d has orthonormal rows. We refer to row(U ) as the planted central subspace. We use the standard notion of subspace coherence, commonly used in matrix completion analysis [11]. Definition 1 (Coherence of the central subspace [11]). Let U ∈ Rr×d satisfy U U ⊤ = Ir . We define the coherence of its row space by µ(U ) :=

d max ∥U:,i ∥22 . r i∈[d]

(5)

Equivalently, with Prow(U ) = U ⊤ U the orthogonal projection onto row(U ), this is the usual subspace coherence µ(row(U )) = dr maxi∈[d] ∥Prow(U ) ei ∥22 with respect to the standard basis (ei )di=1 . Note that small coherence means that no single ambient coordinate carries too much mass from row(U ). Assumption 1 (Finite degree). The link function h is a polynomial of total degree at most ℓ. 4

Throughout, we work with Boolean-hypercube data x ∈ Hd = {±1}d . This setting is natural for sharp KRR analysis since the required low-degree empirical-kernel-matrix approximations are available for uniform spherical/hypercube designs, whereas the corresponding optimal approximation theory for Gaussian designs remains only partially understood and substantially more delicate [29]. On the hypercube, polynomials are naturally reduced modulo the relations x2i = 1. Thus, before taking low-degree truncations, we first replace each polynomial by its multilinear representative. For a polynomial P , let H(P ) denote the unique multilinear polynomial agreeing with P on Hd . Equivalently, on monomials, we have " d # d Y Y αi H xi := xiαi mod 2 , i=1

i=1

and this relation extends linearly to all polynomials on the hypercube. We use the standard convention P that Fq denotes the projection of F onto the degree-q Fourier–Walsh characters, and F≤p := q≤p Fq . Thus H(P )q and H(P )≤p denote the corresponding degree-q component and degree-≤ p truncation of H(P ). See Appendix B.2 for the precise definition of Fourier–Walsh polynomials. On the hypercube, only the values of the ambient polynomial fU∗ (x) = h(U x) on Hd are used. To simplify the notation, we write f ∗ := H(fU∗ ) for its unique multilinear representative. Thus f ∗ (x) = fU∗ (x) for all x ∈ Hd . Throughout the hypercube ∗ analysis, we let fq∗ and f≤p denote the Fourier–Walsh components of this reduced target. The latent link h itself is not reduced and is used below for the Hermite comparison. Next, we define the truncated population AGOP on the hypercube: M≤p =

p X

Mq :=

q=0

p X

(6)

Ex∼τd ∇fq∗ (x)∇fq∗ (x)⊤ .

q=0

Since differentiation lowers Walsh degree by one and distinct Walsh degrees are orthogonal, cross terms between different q’s vanish. Hence we equivalently have ∗ ∗ M≤p = Ex∼τd ∇f≤p (x)∇f≤p (x)⊤ . We will now show that M≤p is closely related to the Gaussian latent AGOP of the degree-≤ p Hermite truncation of h. See Appendix C for the Hermite basis and normalization. Write the Hermite expansion of h under z ∼ N (0, Ir ) as X

h(z) =

α∈Nr : |α|≤ℓ

Letting h≤p (z) :=

aα Heα (z) =

ℓ X

where

hq (z),

q=0

hq (z) :=

aα Heα (z).

α∈Nr : |α|=q

q=0 hq (z), define the degree-≤ p latent gradient covariance

Σp := Ez∼N (0,Ir ) ∇h≤p (z)∇h≤p (z)⊤ ∈ Rr×r .

(7)

The following lemma shows that M≤p is a coherence-controlled perturbation of U ⊤ Σp U : Lemma 1 (Population hypercube-to-Gaussian AGOP comparison). Fix p ∈ {0, . . . , ℓ}. Suppose Assumption 1 holds, and r, ℓ are fixed. Then µ(U ) ∗ 2 M≤p − U ⊤ Σp U op = Od ∥f ∥L2 . (8) d

Learning setup.

We observe training data {(x(i) , yi )}ni=1 with i.i.d.

yi = f ∗ (x(i) ) + εi ,

x(i) ∼ τd ,

i.i.d.

(9)

εi ∼ N (0, σε2 ),

where the noise is independent of the covariates. Let X ∈ Rn×d denote the design matrix with rows x(1) , . . . , x(n) , and let y ∈ Rn denote the response vector. We fit kernel ridge regression with ridge parameter λ ≥ 0 and inner-product kernel ⟨x, x′ ⟩ ′ , (10) K(x, x ) = g d where g : R → R is a kernel profile. The fitted predictor can be written as fˆ(x) = K(x, X) (K(X, X) + λIn )−1 y, where K(x, X) := K(x, x(1) ), . . . , K(x, x(n) ) and K(X, X) := K(x(i) , x(j) ) i,j∈[n] . Throughout the paper, we take the ridge parameter λ = Od (1). We then form the empirical AGOP n X c := 1 M ∇fˆ(x(i) ) ∇fˆ(x(i) )⊤ . n i=1

(11)

Our main results will invoke a combination of the following two regularity assumptions on g. Assumption 2 (Analyticity near 0). There exists ε0 ∈ (0, 1) such that g is analytic on (−ε0 , ε0 ) and satisfies g (k) (0) ≥ 0 for all k ≥ 0. Assumption 3 (Regularity near 1). There exists ε1 ∈ (0, 1) such that g ′ is Lipschitz on (1 − ε1 , 1 + ε1 ). These assumptions cover standard inner-product kernels such as polynomial and exponential kernels, and, on fixed-norm domains, Gaussian kernels.

3.2

Empirical AGOP approximation and subspace recovery

The next theorem gives an operator-norm AGOP approximation with explicit dependence on the subspace coherence µ(U ). We assume that the rank r and largest degree ℓ are fixed as d → ∞. The dependence of the constants on r and ℓ is recorded in the proof that appears in Appendix D. Theorem 2 (AGOP approximation). Suppose Assumptions 1, 2 and 3 hold. Suppose further that P∞ min0≤k≤p g (k) (0) and λ + k=p+1 g (k) (0) are strictly positive. Fix p ∈ {0, . . . , ℓ}, and let n = dp+δ with δ ∈ (0, 1). Then, for every sufficiently small ϵ > 0, the following holds: (12)

c − M≤p ∥op = Od,P Rd (U ) · ∥f ∗ ∥2 + σ 2 , ∥M L2 ε where we set δ

1−δ

1+δ

2+δ

Rd (U ) := d− 2 +ϵ + d− 2 +ϵ + µ(U ) 2 d− 4 +ϵ + µ(U )d− 4 +ϵ + µ(U )2 d− 2 +ϵ .

(13)

In particular, if there exists a constant γ > 0 such that µ(U ) = Od d1/2+δ/4−γ , then, choosing ϵ > 0 c − M≤p ∥op = od,P ∥f ∗ ∥2 + σ 2 . sufficiently small depending on δ and γ yields ∥M ε L2 We next turn the approximation of AGOP into a subspace recovery statement. For row-orthonormal b ∈ Rr×d , we define the usual subspace sine matrix U, U b , U ) := U b (U⊥ )⊤ ∈ Rr×(d−r) . sinΘ (U

b ) and row(U ), and hence are Its singular values are the sines of the principal angles between row(U independent of the choice of U⊥ . We further define sp and ρp which measure the part of M≤p inside and outside row(U ): sp := λmin U M≤p U ⊤ , ρp := M≤p − Prow(U ) M≤p Prow(U ) op , (14) where Prow(U ) = U ⊤ U . The next lemma states the subspace recovery result using sp and ρp . The result follows from the Davis–Kahan sin-theta theorem and the proof appears in Appendix E. c − M≤p Lemma 3 (Eigenspace perturbation). Let εagop := M

b ∈ Rr×d and suppose sp > 0 holds. If U

c, then the estimate holds: has orthonormal rows spanning the top-r eigenspace of M ρp + εagop b, U) sinΘ (U ≤ min 1, 4 . sp op

(15)

Consequently, Lemma 3 reduces subspace recovery to showing (εagop + ρp )/sp = od,P (1). The next corollary verifies this criterion under latent nondegeneracy and weak coherence, by combining Lemma 1, Theorem 2, and Lemma 3. The proof is deferred to Appendix F. Corollary 4 (Subspace recovery under a weak coherence condition). Suppose the assumptions of Theorem 2 hold. Suppose further that there exist constants κ > 0 and γ > 0 such that the following three conditions hold: µ(U ) 2 ∥f ∗ ∥L2 = od (κ). λmin (Σp ) ≥ κ, µ(U ) = Od d1/2+δ/4−γ , d b ∈ Rr×d have orthonormal rows spanning the top-r eigenspace of M c. Then, Let U b, U) sinΘ (U

= od,P κ−1 ∥f ∗ ∥2L2 + σε2

The condition λmin (Σp ) > 0 in Corollary 4 is the AGOP analogue of the nondegeneracy assumptions commonly used in feature-learning analyses of multi-index models [9, 16, 55]: it ensures that the degree-≤ p AGOP contains all relevant directions. The coherence condition is a sufficient technical condition. It is substantially weaker than the bounded or polylogarithmic incoherence assumptions common in matrix completion analysis [11]. The experiments in Section 5 suggest that AGOP-based subspace recovery succeeds under broader conditions than those guaranteed by Corollary 4.

Implications for feature learning

Recursive Feature Machines (RFM) are iterative kernel methods proposed to enable feature learning in kernel machines: they alternate between fitting a kernel predictor and updating the kernel metric using the empirical AGOP of that predictor [46, 47]. This section shows how one RFM step turns AGOP recovery into feature learning. Although the first KRR fit from the isotropic metric M1 = Id mainly captures the low-degree component of the target, its empirical AGOP already recovers the central subspace row(U ). RFM then uses this AGOP to update the metric, producing an anisotropic rescaling: relevant directions are amplified, while orthogonal directions remain essentially unchanged. We consider inner-product kernels of the form ⊤ ⟨M 1/2 x, M 1/2 x′ ⟩ x M x′ ′ KM (x, x ) = g =g , (16) d d where M ∈ Rd×d is positive semidefinite. Therefore, KRR with kernel KM can be interpreted as ordinary inner-product-kernel KRR applied to the transformed inputs M 1/2 x. Starting from M1 = Id , and fixing a safeguard parameter η > 0 and ridge parameter λ > 0, RFM alternates between 7

• Step 1 (KRR): fˆt (x) = KMt (x, X) (KMt (X, X) + λIn )−1 y. • Step 2 (metric update): Mt+1 =

d ct + ηId ) tr(M

ct + ηId ), (M

n X ct := 1 M ∇fˆt (x(i) ) ∇fˆt (x(i) )⊤ . n i=1

The trace normalization fixes the metric scale across iterations, while the safeguard ηId , as in [56], prevents small-AGOP directions from being driven to zero. At the first iteration, fˆ1 is exactly the standard KRR predictor. The next theorem, adapted from Theorem 4 of [37], records that in the polynomial-sample regime it recovers only the low-degree part of the target. The proof is deferred to Appendix G. Theorem 5 (Adapted from Theorem 4 of [37]). Under the assumptions of Theorem 2, the first-step KRR predictor fˆ1 satisfies ∗ ∥fˆ1 − f≤p ∥2L2 = od,P (1) ∥f ∗ ∥2L2 + σε2 . (17) Combining Theorem 5 with Corollary 4 yields a prediction–representation separation. Indeed, if the target function f ∗ has degree p∗ with p∗ > p, then in the regime n = dp+δ , the empirical AGOP recovers the central subspace, while the prediction error of the first KRR fit can remain large: 2 ∗ ∥fˆ1 − f ∗ ∥2L2 = f>p + od,P (1) ∥f ∗ ∥2L2 + σε2 . L 2

(18)

In the remainder of this section, suppose that all assumptions of Corollary 4 hold. Then, more c1 already identifies the central subspace. importantly for our purposes, Corollary 4 shows that the AGOP M The next proposition shows that the first RFM metric update turns this recovered geometry into an anisotropic rescaling of the inputs. The proof is deferred to Appendix H. Proposition 6. Set the safeguard parameter η = dζ · max{d−δ/2 , d−(1−δ)/2 } for a small constant ζ > 0. For an independent sample x ∼ τd , define p √ x b := U ⊤ Σp + ηIr U x + η (U⊥ )⊤ U⊥ x. (19) c1 + ηId ). Then the following estimate holds Let cη := tr(M 1/2

M2 x −

d/cη x b L = od,P 2

q d/cη ∥b x∥L2 ,

(20)

where cη = (1 + od,P (1)) ηd. Writing z0 := U x and z1 := U⊥ x, Proposition 6 indicates q 1/2 ⊤ ⊤ −1 M2 x = (1 + od,P (1)) U η Σp + Ir z0 + (U⊥ ) z1 .

(21)

Thus the second p iterate applies an anisotropic linear map within the relevant subspace, with singular values at least η −1 Σp + Ir , while leaving the orthogonal coordinates at unit scale. In particular, since p Σp ⪰ κIr then every relevant direction is amplified by at least 1 + κ/η. This provides a mechanism by which one RFM update performs feature learning.

Conjectural prediction for the second fit. The representation above suggests a heuristic analogy with the spiked-covariate model of [23]. Motivated by this analogy, we conjecture that the relevant sample √ complexity for the second KRR fit is governed by the effective dimension deff := d η. Equivalently, one log d expects the polynomial threshold to shift from p to peff = p · log(d ≥ p. In particular, this heuristic eff ) ˆ predicts that the second-step KRR predictor f2 with independent draws of input data, should satisfy ∗ ∥fˆ2 − f≤⌊p ∥2 = od,P (1) ∥f ∗ ∥2L2 + σε2 . eff ⌋ L2

(22)

We leave this as a conjectural extension. The rigorous contribution of this section is the anisotropic rescaling result above, which explains how one RFM update can convert recovered subspace information into a more favorable geometry for the next KRR fit.

Numerical results

We evaluate two phenomena highlighted by our analysis: recovery of the central subspace from the AGOP, and the empirical improvement produced by iterating RFM. Target functions.

Let z = U x with U ∈ Rr×d . We consider the two link functions

1 h(z) = He1 (z1 ) + √ He4 (z1 ) (L1), 24

h(z) = He1 (z1 )He1 (z2 ) +

1 He2 (z1 )He2 (z2 ) (L2). 2

The first is a single-index model (r = 1), and the second is an index-two model (r = 2). In both cases, the low-degree component already contains the full subspace information, while the higher-degree term increases prediction difficulty. Input data.

We consider two input distributions: x ∼ Unif({−1, 1}d ),

x ∼ N (0, Id ).

The Boolean distribution matches the design studied in our theoretical analysis. For both input distributions, U is formed by taking the first r rows of a Haar-distributed orthogonal matrix. We also consider a Boolean design, x ∼ Unif({−1, 1}d ), with a sparse matrix U . This setting does not satisfy the incoherence condition in Corollary 4 hence it is for testing how the method behaves when this condition fails. We construct U so that its rows are sparse and have pairwise disjoint supports. On each support, we draw independent random entries and then normalize the resulting vector to have unit Euclidean norm. The support size is set to d3/10 ≈ 4 when d = 100. Experimental protocol.

We fix the ambient dimension to d = 100 and generate the labels by ε ∼ N (0, 0.01).

y = h(U x) + ε,

The test loss, specifically mean squared error, is evaluated on 5, 000 test points. For each sample-size √ exponent α, we set n = ⌊dα ⌋. For the Gaussian kernel we use bandwidth d, for the Laplace kernel d, and the ridge parameter is λ = 10−6 . The RFM safeguard is η = 10−2 · d. We run RFM with Gaussian and Laplace kernels for five iterations and report averages over 10 independent trials. We report the experimental results for the Gaussian kernel with hypercube data and dense U in Figs. 1 and 2, corresponding to the target functions (L1) and (L2), respectively. The remaining plots are provided in Appendix A. The results show that, for KRR, corresponding to iteration 0, increasing the sample-size exponent leads to improved subspace recovery: the principal angle decreases and eventually becomes small, even though the test loss remains relatively large; see the first and second panels. After additional iterations of RFM, both the test loss and the principal angle improve substantially. These findings support our 9

theoretical results on subspace recovery and are consistent with the feature-learning mechanism suggested by our RFM analysis.

Sample-size exponent

1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.1 1.0 0.9 0.8

||sin (U, U)||op

Test loss

1.0

1.5

Iteration

Third eigenvalue

10 2

0.4

10 3

0.2 0

Iteration

0.0

100 10 1

0.6

0.5 1

Second eigenvalue

0.8

1.0

First eigenvalue

Iteration

Figure 1: For target function (L1), we train RFM with a Gaussian kernel on hypercube data for five iterations. From left to right, the panels display the test loss, the sine of the largest principal angle, and the largest, second-largest and third-largest eigenvalues of the AGOP.

Sample-size exponent

2.2 2.1 2.0 1.9 1.8 1.7 1.6 1.5

||sin (U, U)||op

Test loss

1.0

10 1 10 2

0.4

0.5

Iteration

Third eigenvalue

0.6

1.0

Second eigenvalue

0.8

1.5

First eigenvalue

10 3

0.2 0

Iteration

0.0

Iteration

Figure 2: For target function (L2), we train RFM with a Gaussian kernel on hypercube data for five iterations. From left to right, the panels display the test loss, the sine of the largest principal angle, and the largest, second-largest and third-largest eigenvalues of the AGOP.

Conclusion

We studied whether a fitted predictor can reveal low-dimensional structure before it predicts accurately. For polynomial multi-index models on the Boolean hypercube, we showed that the empirical AGOP of a KRR predictor concentrates around the population AGOP of a low-degree truncation of the target. If this population AGOP is nondegenerate on the planted subspace and the subspace satisfies a weak coherence condition, then the leading eigenspace of the empirical AGOP recovers the planted central subspace. This gives a separation between prediction and representation. In the regime n = dp+δ , KRR may learn only the degree-≤ p part of the target, so its prediction error can remain large when higher-degree components are present. Nevertheless, the AGOP of the fitted predictor can already identify the central subspace once the degree-≤ p population AGOP is nondegenerate on that subspace. Thus, representation recovery is governed by the first degree at which this nondegeneracy holds, rather than by the highest degree needed for accurate prediction. Finally, we showed that one RFM update uses the recovered AGOP to rescale the input anisotropically: directions in the central subspace are amplified, while orthogonal directions remain essentially unchanged. This provides a mechanism by which iterative AGOP-based kernel methods can turn early representation recovery into improved subsequent fits in multi-index models.

References [1] Emmanuel Abbe, Enric Boix Adsera, and Theodor Misiakiewicz. “Sgd learning on neural networks: leap complexity and saddle-to-saddle dynamics”. In: The Thirty Sixth Annual Conference on Learning Theory. PMLR. 2023, pp. 2552–2623. [2] Emmanuel Abbe, Enric Boix Adsera, and Theodor Misiakiewicz. “The merged-staircase property: a necessary and nearly sufficient condition for sgd learning of sparse functions on two-layer neural networks”. In: Conference on Learning Theory. PMLR. 2022, pp. 4782–4887. [3] George E Andrews and Richard Askey. “Classical orthogonal polynomials”. In: Polynômes Orthogonaux et Applications: Proceedings of the Laguerre Symposium held at Bar-le-Duc, October 15–18, 1984. Springer. 2006, pp. 36–62. [4] Luca Arnaboldi, Yatin Dandi, Florent Krzakala, Bruno Loureiro, Luca Pesce, and Ludovic Stephan. “Online learning and information exponents: On the importance of batch size, and time/complexity tradeoffs”. In: arXiv preprint arXiv:2406.02157 (2024). [5] Luca Arnaboldi, Yatin Dandi, Florent Krzakala, Luca Pesce, and Ludovic Stephan. “Repetita iuvant: Data repetition allows sgd to learn high-dimensional multi-index functions”. In: arXiv preprint arXiv:2405.15459 (2024). [6] Gerard Ben Arous, Reza Gheissari, and Aukosh Jagannath. “Online stochastic gradient descent on non-convex losses from high-dimensional inference”. In: Journal of Machine Learning Research 22.106 (2021), pp. 1–51. [7] Jimmy Ba, Murat A Erdogdu, Taiji Suzuki, Zhichao Wang, Denny Wu, and Greg Yang. “Highdimensional asymptotics of feature learning: How one gradient step improves the representation”. In: Advances in Neural Information Processing Systems 35 (2022), pp. 37932–37946. [8] Daniel Beaglehole, David Holzmüller, Adityanarayanan Radhakrishnan, and Mikhail Belkin. “xRFM: Accurate, scalable, and interpretable feature learning models for tabular data”. In: The Fourteenth International Conference on Learning Representations. 2026. url: https://openreview.net/ forum?id=wHuVdpnUFp. [9]

Alberto Bietti, Joan Bruna, and Loucas Pillaud-Vivien. “On learning gaussian multi-index models with gradient flow”. In: arXiv preprint arXiv:2310.19793 (2023).

[10]

Alberto Bietti, Joan Bruna, Clayton Sanford, and Min Jae Song. “Learning single-index models with shallow neural networks”. In: Advances in Neural Information Processing Systems 35 (2022), pp. 9768–9783.

[11]

Emmanuel Candes and Benjamin Recht. “Exact matrix completion via convex optimization”. In: Communications of the ACM 55.6 (2012), pp. 111–119.

[12]

Theodore S Chihara. An introduction to orthogonal polynomials. Courier Corporation, 2011.

[13]

R Dennis Cook. “Fisher lecture: Dimension reduction in regression”. In: (2007).

[14]

Alex Damian, Jason D Lee, and Joan Bruna. “The Generative Leap: Tight Sample Complexity for Efficiently Learning Gaussian Multi-Index Models”. In: The Thirty-ninth Annual Conference on Neural Information Processing Systems.

[15]

Alex Damian, Eshaan Nichani, Rong Ge, and Jason D. Lee. “Smoothing the Landscape Boosts the Signal for SGD: Optimal Sample Complexity for Learning Single Index Models”. In: Thirty-seventh Conference on Neural Information Processing Systems. 2023. url: https://openreview.net/ forum?id=73XPopmbXH.

[16]

Alexandru Damian, Jason Lee, and Mahdi Soltanolkotabi. “Neural networks can learn representations with gradient descent”. In: Conference on Learning Theory. PMLR. 2022, pp. 5413–5452.

[17]

Yatin Dandi, Florent Krzakala, Bruno Loureiro, Luca Pesce, and Ludovic Stephan. “How two-layer neural networks learn, one (giant) step at a time”. In: Journal of Machine Learning Research 25.349 (2024), pp. 1–65. 11

[18]

Yatin Dandi, Emanuele Troiani, Luca Arnaboldi, Luca Pesce, Lenka Zdeborova, and Florent Krzakala. “The Benefits of Reusing Batches for Gradient Descent in Two-Layer Networks: Breaking the Curse of Information and Leap Exponents”. In: International Conference on Machine Learning. PMLR. 2024, pp. 9991–10016.

[19]

Noureddine El Karoui. “The spectrum of kernel random matrices”. In: The Annals of Statistics 38.1 (2010), pp. 1–50. doi: 10.1214/09-AOS715.

[20]

Kenji Fukumizu, Francis R Bach, and Michael I Jordan. “Kernel dimension reduction in regression”. In: The Annals of Statistics 37.4 (2009), pp. 1871–1905. doi: 10.1214/08-AOS637.

[21]

Kenji Fukumizu and Chenlei Leng. “Gradient-based kernel dimension reduction for regression”. In: Journal of the American Statistical Association 109.505 (2014), pp. 359–370.

[22]

Behrooz Ghorbani, Song Mei, Theodor Misiakiewicz, and Andrea Montanari. “Linearized two-layers neural networks in high dimension”. In: The Annals of Statistics 49.2 (2021).

[23]

Behrooz Ghorbani, Song Mei, Theodor Misiakiewicz, and Andrea Montanari. “When do neural networks outperform kernel methods?” In: Advances in Neural Information Processing Systems 33 (2020), pp. 14820–14830.

[24]

Wolfgang Härdle and Thomas M Stoker. “Investigating smooth multiple regression by the method of average derivatives”. In: Journal of the American statistical Association 84.408 (1989), pp. 986–995.

[25]

Marian Hristache, Anatoli Juditsky, Jörg Polzehl, and Vladimir Spokoiny. “Structure adaptive approach for dimension reduction”. In: Annals of Statistics (2001), pp. 1537–1566.

[26]

Shuo Huang, Hippolyte Labarrière, Ernesto De Vito, Tomaso Poggio, and Lorenzo Rosasco. “Learning Multi-Index Models with Hyper-Kernel Ridge Regression”. In: arXiv preprint arXiv:2510.02532 (2025).

[27]

Hidehiko Ichimura. “Semiparametric least squares (SLS) and weighted SLS estimation of single-index models”. In: Journal of econometrics 58.1-2 (1993), pp. 71–120.

[28]

Arthur Jacot, Franck Gabriel, and Clément Hongler. “Neural tangent kernel: Convergence and generalization in neural networks”. In: Advances in neural information processing systems. 2018, pp. 8571–8580.

[29]

Chiraag Kaushik, Justin Romberg, and Vidya Muthukumar. “A general technique for approximating high-dimensional empirical kernel matrices”. In: arXiv preprint arXiv:2511.03892 (2025).

[30]

Alex Kokot, Anand Hemmady, Vydhourie Thiyageswaran, and Marina Meila. “Local EGOP for Continuous Index Learning”. In: arXiv preprint arXiv:2601.07061 (2026).

[31]

Michel Ledoux and Michel Talagrand. Probability in Banach Spaces: isoperimetry and processes. Springer Science & Business Media, 2013.

[32]

Jason D Lee, Kazusato Oko, Taiji Suzuki, and Denny Wu. “Neural network learns low-dimensional polynomials with sgd near the information-theoretic limit”. In: Advances in Neural Information Processing Systems 37 (2024), pp. 58716–58756.

[33]

Ker-Chau Li. “On principal Hessian directions for data visualization and dimension reduction: another application of Stein’s lemma”. In: Journal of the American Statistical Association 87.420 (1992), pp. 1025–1039.

[34]

Ker-Chau Li. “Sliced Inverse Regression for Dimension Reduction”. In: Journal of the American Statistical Association 86.414 (1991), pp. 316–327. issn: 01621459, 1537274X.

[35]

Qian Lin, Zhigen Zhao, and Jun S Liu. “ON CONSISTENCY AND SPARSITY FOR SLICED INVERSE REGRESSION IN HIGH DIMENSIONS”. In: The Annals of Statistics 46.2 (2018), pp. 580–610.

[36]

Yanyuan Ma and Liping Zhu. “A review on dimension reduction”. In: International Statistical Review 81.1 (2013), pp. 134–150.

[37]

Song Mei, Theodor Misiakiewicz, and Andrea Montanari. “Generalization error of random feature and kernel methods: hypercontractivity and kernel matrix concentration”. In: Applied and Computational Harmonic Analysis 59 (2022), pp. 3–84.

[38]

Behrad Moniri, Donghwan Lee, Hamed Hassani, and Edgar Dobriban. “A theory of non-linear feature learning with one gradient step in two-layer neural networks”. In: Proceedings of the 41st International Conference on Machine Learning. 2024, pp. 36106–36159.

[39]

Andrea Montanari and Zihao Wang. “Phase transitions for feature learning in neural networks”. In: arXiv preprint arXiv:2602.01434 (2026).

[40]

Alireza Mousavi-Hosseini, Sejun Park, Manuela Girotti, Ioannis Mitliagkas, and Murat A Erdogdu. “Neural Networks Efficiently Learn Low-Dimensional Representations with SGD”. In: The Eleventh International Conference on Learning Representations. 2023.

[41]

[42]

Alireza Mousavi-Hosseini, Denny Wu, and Murat A Erdogdu. “Learning Multi-Index Models with Neural Networks via Mean-Field Langevin Dynamics”. In: The Thirteenth International Conference on Learning Representations. 2025. url: https://openreview.net/forum?id=WHhZv8X5zF.

[43]

Ryan O’Donnell. Analysis of boolean functions. Cambridge University Press, 2014.

[44]

Parthe Pandit, Zhichao Wang, and Yizhe Zhu. “Universality of kernel random matrices and kernel regression in the quadratic regime”. In: Journal of Machine Learning Research 26.224 (2025), pp. 1– 73.

[45]

James L Powell, James H Stock, and Thomas M Stoker. “Semiparametric estimation of index coefficients”. In: Econometrica: Journal of the Econometric Society (1989), pp. 1403–1430.

[46]

Adityanarayanan Radhakrishnan, Daniel Beaglehole, Parthe Pandit, and Mikhail Belkin. “Mechanism for feature learning in neural networks and backpropagation-free machine learning models”. In: Science 383.6690 (2024), pp. 1461–1467. doi: 10.1126/science.adi5639. eprint: https://www. science.org/doi/pdf/10.1126/science.adi5639. url: https://www.science.org/doi/abs/ 10.1126/science.adi5639.

[47]

Adityanarayanan Radhakrishnan, Mikhail Belkin, and Dmitriy Drusvyatskiy. “Linear recursive feature machines provably recover low-rank matrices”. In: Proceedings of the National Academy of Sciences 122.13 (2025), e2411325122.

[48]

Feng Ruan, Keli Liu, and Michael Jordan. “A Compositional Kernel Model for Feature Learning”. In: arXiv preprint arXiv:2509.14158 (2025).

[49]

Berfin Simsek, Amire Bendjeddou, and Daniel Hsu. “Learning Gaussian Multi-Index Models with Gradient Flow: Time Complexity and Directional Convergence”. In: The 28th International Conference on Artificial Intelligence and Statistics. 2025. url: https://openreview.net/forum?id=wYfOdkKnsr.

[50]

Gabor Szeg. Orthogonal polynomials. Vol. 23. American Mathematical Soc., 1939.

[51]

Konstantinos Christopher Tsiolis, Alireza Mousavi-Hosseini, and Murat A Erdogdu. “From Information to Generative Exponent: Learning Rate Induces Phase Transitions in SGD”. In: The Thirty-ninth Annual Conference on Neural Information Processing Systems (2025).

[52]

Yingcun Xia. “A Constructive Approach to the Estimation of Dimension Reduction Directions”. In: The Annals of Statistics (2007), pp. 2654–2690.

[53]

Yingcun Xia, Howell Tong, Wai Keung Li, and Li-Xing Zhu. “An adaptive estimation of dimension reduction space”. In: Journal of the Royal Statistical Society Series B: Statistical Methodology 64.3 (2002), pp. 363–410.

[54]

Gan Yuan, Mingyue Xu, Samory Kpotufe, and Daniel Hsu. “Efficient estimation of the central mean subspace via smoothed gradient outer products”. In: SIAM Journal on Mathematics of Data Science 7.3 (2025), pp. 1241–1264.

[55]

Bohan Zhang, Zihao Wang, Hengyu Fu, and Jason D. Lee. “Neural Networks Learn Generic MultiIndex Models Near Information-Theoretic Limit”. In: The Fourteenth International Conference on Learning Representations. 2026. url: https://openreview.net/forum?id=2Q0U2rV2Jz.

[56]

Libin Zhu, Damek Davis, Dmitriy Drusvyatskiy, and Maryam Fazel. “Iteratively reweighted kernel machines efficiently learn sparse functions”. In: arXiv preprint arXiv:2505.08277 (2025).

[57]

Libin Zhu, Damek Davis, Dmitriy Drusvyatskiy, and Maryam Fazel. “Spectral norm bound for the product of random Fourier-Walsh matrices”. In: arXiv preprint arXiv:2504.03148 (2025).

Appendix Contents A Additional numerical results

B Preliminaries B.1 Preliminaries on orthogonal polynomials . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.2 Fourier expansion on the boolean Hypercube . . . . . . . . . . . . . . . . . . . . . . . . . B.3 Isometric properties of orthogonal polynomials . . . . . . . . . . . . . . . . . . . . . . . . B.4 Approximating kernels by orthogonal polynomials . . . . . . . . . . . . . . . . . . . . . . .

19 19 20 21 22

C Proof of Lemma 1

D A General AGOP Approximation Result and Proof of Theorem 2 D.1 Proof of Theorem 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31 34

E Proof of Lemma 3

F Proof of Corollary 4

G Proof of Theorem 5

H Proof of Proposition 6

Proof for preliminaries I.1 Proof of Lemma 13 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I.2 Proof of Lemma 14 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I.3 Proof of Lemma 16 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I.4 Proof of Theorem 17 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I.5 Proof of Proposition 18 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I.6 Proof of Proposition 19 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I.7 Proof of Lemma 20 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I.8 Proof of Proposition 21 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I.9 Proof of Proposition 22 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I.10 Proof of Theorem 23 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42 42 43 45 46 47 48 51 51 53 54

J Proof of Theorem 25 J.1 Proof of Item (a) of Claim 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . J.2 Proof of Item (b) of Claim 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

55 56 57

K Proof of lemmas and propositions for the main results K.1 Proof of Lemma 27 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K.2 Proof of Lemma 28 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K.3 Proof of Proposition 29 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K.4 Proof of Proposition 30 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K.5 Proof of Claim 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K.6 Proof of Proposition 31 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K.7 Proof of Proposition 32 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K.8 Proof of Proposition 33 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K.9 Proof of Proposition 34 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K.10 Proof of Proposition 35 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K.11 Proof of Proposition 36 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K.12 Proof of Claim 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

60 60 63 64 65 68 68 70 72 74 75 77 78

K.13 Proof of Proposition 37 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K.14 Proof of Proposition 38 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K.15 Proof of Claim 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

79 81 82

L Product of random Fourier-Walsh matrices L.1 Proof of Proposition 40 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . L.2 Proof of Proposition 41 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

83 90 91

M Technical lemmas and propositions

Additional numerical results

In this section, we report the complementary experimental results for Section 5. The experiments follow the same protocol as in the main text. In each plot, iteration 0 corresponds to KRR, while the later iterations correspond to RFM updates. Across the additional settings, we observe the same qualitative behavior: increasing the sample-size exponent improves the initial subspace estimate, and subsequent RFM iterations further reduce both the test loss and the sine of the largest principal angle. The evolution of the leading AGOP eigenvalues is also consistent with the emergence of the relevant low-dimensional structure. Gaussian kernel.

Sample-size exponent

1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.1 1.0 0.9 0.8

||sin (U, U)||op

Test loss

1.0

10 2 10 3

0.2 0

Iteration

0.0

100 10 1

0.4

0.5

Iteration

Third eigenvalue

0.6

1.0

Second eigenvalue

0.8

1.5

First eigenvalue

Iteration

Figure 3: For target function (L1), we train RFM with a Gaussian kernel on dense Gaussian data for five iterations. From left to right, the panels display the test loss, the sine of the largest principal angle, and the largest, second-largest and third-largest eigenvalues of the AGOP.

Sample-size exponent

2.2 2.1 2.0 1.9 1.8 1.7 1.6 1.5

Test loss

||sin (U, U)||op

2.0

1.0

10 1 10 2

0.4

0.5

Iteration

Third eigenvalue

0.6

1.0 1

Second eigenvalue

0.8

1.5

First eigenvalue

0.2 0

Iteration

0.0

10 3 0

Iteration

Figure 4: For target function (L2), we train RFM with a Gaussian kernel on dense Gaussian data for five iterations. From left to right, the panels display the test loss, the sine of the largest principal angle, and the largest, second-largest and third-largest eigenvalues of the AGOP.

Sample-size exponent

1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.1 1.0 0.9 0.8

Test loss

Iteration

||sin (U, U)||op

1.25

1.0

1.00

0.8

0.75

0.6

0.50

0.4

0.25

0.2

Iteration

0.0

First eigenvalue

Second eigenvalue

Third eigenvalue 10 1 10 2 10 3

Iteration

Figure 5: For target function (L1), we train RFM with a Gaussian kernel on sparse hypercube data for five iterations. From left to right, the panels display the test loss, the sine of the largest principal angle, and the largest, second-largest and third-largest eigenvalues of the AGOP.

Sample-size exponent

2.2 2.1 2.0 1.9 1.8 1.7 1.6 1.5

||sin (U, U)||op

Test loss

Iteration

1.0

1.00

0.8

0.75

0.6

0.50

0.4

0.25

0.2

Iteration

0.0

First eigenvalue

Second eigenvalue

Third eigenvalue 10 1 10 2 10 3

Iteration

Figure 6: For target function (L2), we train RFM with a Gaussian kernel on sparse hypercube data for five iterations. From left to right, the panels display the test loss, the sine of the largest principal angle, and the largest, second-largest and third-largest eigenvalues of the AGOP. Laplacian kernel.

Sample-size exponent

1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.1 1.0 0.9 0.8

Test loss

||sin (U, U)||op

2.0

1.0

0.5 2

Iteration

Third eigenvalue 100 10 1

0.6

1.0

Second eigenvalue

0.8

1.5

First eigenvalue

Iteration

0.4

10 2

0.2

10 3

0.0

Iteration

Figure 7: For target function (L1), we train RFM with a Laplacian kernel on dense Gaussian data for five iterations. From left to right, the panels display the test loss, the sine of the largest principal angle, and the largest, second-largest and third-largest eigenvalues of the AGOP.

Sample-size exponent

2.2 2.1 2.0 1.9 1.8 1.7 1.6 1.5

Test loss

||sin (U, U)||op

2.0

1.0

0.5 2

Iteration

Third eigenvalue

Iteration

100 10 1

0.6

1.0 1

Second eigenvalue

0.8

1.5

First eigenvalue

0.4

10 2

0.2

10 3

0.0

Iteration

Figure 8: For target function (L2), we train RFM with a Laplacian kernel on dense Gaussian data for five iterations. From left to right, the panels display the test loss, the sine of the largest principal angle, and the largest, second-largest and third-largest eigenvalues of the AGOP.

Sample-size exponent

1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.1 1.0 0.9 0.8

||sin (U, U)||op

Test loss

1.0

1.5

Iteration

Third eigenvalue 10 1

0.6 0.4

0.5 1

Second eigenvalue

0.8

1.0

First eigenvalue

10 3

0.2 0

Iteration

0.0

Iteration

Figure 9: For target function (L1), we train RFM with a Laplacian kernel on dense hypercube data for five iterations. From left to right, the panels display the test loss, the sine of the largest principal angle, and the largest, second-largest and third-largest eigenvalues of the AGOP.

Sample-size exponent

2.2 2.1 2.0 1.9 1.8 1.7 1.6 1.5

Test loss

||sin (U, U)||op

2.0

1.0

0.5 2

Iteration

Third eigenvalue

Iteration

100 10 1

0.6

1.0

Second eigenvalue

0.8

1.5

First eigenvalue

0.4

10 2

0.2

10 3

0.0

Iteration

Figure 10: For target function (L2), we train RFM with a Laplacian kernel on dense hypercube data for five iterations. From left to right, the panels display the test loss, the sine of the largest principal angle, and the largest, second-largest and third-largest eigenvalues of the AGOP.

Sample-size exponent

1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.1 1.0 0.9 0.8

||sin (U, U)||op

Test loss

Iteration

1.0

0.8

0.5

0.4

First eigenvalue

Second eigenvalue

Third eigenvalue 10 1

0.6

10 3

0.2 0

Iteration

0.0

Iteration

Figure 11: For target function (L1), we train RFM with a Laplacian kernel on sparse hypercube data for five iterations. From left to right, the panels display the test loss, the sine of the largest principal angle, and the largest, second-largest and third-largest eigenvalues of the AGOP.

Sample-size exponent

2.2 2.1 2.0 1.9 1.8 1.7 1.6 1.5

||sin (U, U)||op

Test loss

Iteration

1.0

0.8

0.5

0.4

First eigenvalue

Second eigenvalue

Third eigenvalue 10 1

0.6

10 3

0.2 0

Iteration

0.0

Iteration

Figure 12: For target function (L2), we train RFM with a Laplacian kernel on sparse hypercube data for five iterations. From left to right, the panels display the test loss, the sine of the largest principal angle, and the largest, second-largest and third-largest eigenvalues of the AGOP.

Preliminaries

B.1

Preliminaries on orthogonal polynomials

In this section, we record notation and preliminaries on orthogonal polynomials, following the standard monographs on the subject [3, 12, 50]. Consider a measure space (D, A, µ), where D is a set, A is a σ-algebra, and µ is a measure on (D, A). We let L2 (D, µ) denote the Hilbert space of µ square-integrable functions on D equipped with the usual inner product and the induced norm Z p ⟨f, g⟩ = f g dµ ∥f ∥ = ⟨f, f ⟩. We will suppress the symbol D from L2 (D, µ) when it is clear from context. Every Hilbert space L2 (µ) admits an orthonormal basis {ϕj }j∈J , meaning ϕj are unit norm, pairwise orthogonal, and their linear span is dense in L2 (µ). A favorable situation occurs when J is countable, in which case L2 (µ) is called separable. Any function f ∈ L2 (µ) in a separable Hilbert space can be expanded in the orthogonal basis: f−

m X ⟨f, ϕi ⟩ϕi → 0

as m ↗ |J|.

i=0

Here, we identify J with the contiguous subset of natural numbers N starting at zero. We will primarily focus on the following two examples of separable Hilbert spaces along with orthonormal bases.

B.2

Fourier expansion on the boolean Hypercube

Let H = {−1, 1}d be the hypercube equipped with the uniform measure τd . Then the inner product between any two functions f, g : Hd → R is simply d

⟨f, g⟩ =

1 X f (x)g(x). 2d d x∈H

An orthonormal basis on L2 (Hd ) is furnished by the multi-linear monomials (called Fourier-Walsh) α

x :=

d Y

for each α ∈ {0, 1}d .

i xα i

i=1

Pd We will denote the degree of the monomial xα by the symbol |α| := i=1 αi . Notice that there is a one-to-one correspondence between binary vectors α ∈ {0, 1}d and subsets S ⊂ [d] (their support). We will therefore often abuse notation and treat α as both a vector and a set whenever convenient. Thus any polynomial f on the hypercube Hd can be expanded in the monomial basis: X f (x) = bα xα with bα = ⟨f, xα ⟩, α∈{0,1}d

where bα is called the Fourier coefficient of f indexed by α. The p-truncation of f is then defined to be the truncated series X f≤p (x) = bα x α . (23) α∈{0,1}d : |α|≤p

Equivalently, f≤p is the projection of f in L2 (Hd ) onto the span of all Fourier-Walsh monomials xα with |α| ≤ p. The functions fp and f>p are defined in an obvious way. The Fourier coefficients bα have a convenient interpretation in terms of discrete derivative of f . Namely, the discrete derivative of any function g : Hd → R in direction xi is defined to be Di g(x) =

g(x) − g(x1 , . . . , xi−1 , −xi , xi+1 , . . . , xd ) . 2xi

(24)

The discrete derivative with respect to a set S = {i1 , . . . , ik } ⊂ [d] is then defined by iterating: DS f (x) = Di1 Di2 . . . Dik f (x). If G := H(g) : Rd → R denotes the multilinear extension of g, then this normalization agrees with the ambient derivatives on the hypercube: x ∈ Hd , i ∈ [d].

Di g(x) = ∂i G(x), More generally, for every S ⊂ [d], DS g(x) = ∂S G(x),

x ∈ Hd ,

Q where ∂S := i∈S ∂i . Interestingly, Fourier coefficients correspond precisely to expectations of discrete derivatives: bS = E[DS f ].

(25)

See for example [43, Section 2.2] for details. Therefore, Fourier coefficients measure the sensitivity of f to coordinate perturbations. Consider now Rd equipped with some probability measure µd . We will often have to control the Lq -norm ∥f ∥Lq (µd ) = [E|f |q ]1/q of a polynomial f on Rd . In general, one would expect that the ratio 20

between the Lq (µd ) and L2 (µd ) norms depends strongly on the dimension d. Interestingly, for a broad class of measures µ and functions f this is not the case. Indeed, the so-called hypercontractivity inequality ensures that any polynomial f : Rd → R of degree at most ℓ satisfies the inequality ∥f ∥Lq (µd ) ≤ (q − 1)ℓ/2 ∥f ∥L2 (µd )

∀q ≥ 2,

(26)

where µd can be the standard Gaussian measure [31, Chapter 3.2] on Rd or the uniform measure on the hypercube Hd [43, Chapter 9]. The inequality (26) is widely used in probability and theoretical computer science, and we will use it heavily here as well.

B.3

Isometric properties of orthogonal polynomials

In this subsection, we restate several isometric properties of covariance-like matrices induced by orthogonal polynomials. These results are established in [56]; we include them here for convenience, with minor notational changes, and omit the proofs. Fix a probability space (D, A, µd ) with D ⊂ Rd and a set {ϕj }j∈S̄ of orthonormal polynomials with respect to µd , indexed by some set S̄. For any finite set S ⊂ S̄ and a point x ∈ D, we define the concatenated vector ϕS (x) := (ϕj (x))j∈S ∈ R|S| . Given a sequence of points X = (x(1) , . . . , x(n) ), we stack ϕS (x(i) ) as rows to form the matrix ΦS (X) = [ϕj (x(i) )] ∈ Rn×|S| . The rows of ΦS (X) are indexed by the data points and the columns by the polynomials in S. To simplify notation, we will often omit X from ΦS (X) when it is clear from context. Throughout, we assume that the data points x(1) , . . . , x(n) are sampled independently from µd . We impose the following hypercontractivity assumption throughout this subsection, paralleling (26) in our two running examples. More precisely, we assume hypercontractivity of the product measure µd × µd on Rd × Rd . Hypercontractivity of the original measure µd then follows immediately. Assumption 4 (Hypercontractivity). There exist constants Cℓ,q > 0, indexed by integers ℓ, q ∈ N, such that any polynomial f on Rd × Rd of degree at most ℓ satisfies ∥f ∥Lq (µd ×µd ) ≤ Cℓ,q · ∥f ∥L2 (µd ×µd )

∀q ≥ 2.

In particular, Assumption 4 holds in our running example—uniform on the hypercube—with Cℓ,q = (q − 1)ℓ/2 , since the product uniform measure uniform on the hypercube is again uniform on the hypercube. We now record the relevant isometry properties of the matrix ΦS in the high-dimensional regime n, d → ∞. Lemma 7 (Norm control). Suppose Assumption 4 holds. Consider the regime n = dp+δ where δ ∈ (0, 1) is a constant. Fix a set S ⊆ S̄ indexing polynomials of degree at most ℓ. Then the following hold: 1. If |S| ≤ Cdp+δ0 with δ0 ∈ (−∞, δ), then √ ∥ΦS ∥op = Od,P ( n). 2. If |S| ≥ Cdp+δ0 with δ0 ∈ (δ, ∞), then p ∥ΦS ∥op = Od,P ( |S|). The asymptotic behavior of ΦS depends on the relative scale of n and |S|, namely n ≍ dq with q>0 or q<0. | {z } | {z } |S| Case I

Case II

The following two theorems show, respectively, that in Case I the empirical covariance matrix n1 Φ⊤ S ΦS is 1 asymptotically equal to the identity I|S| , while in Case II the Gram matrix |S| ΦS Φ⊤ is asymptotically S equal to the identity In . 21

Theorem 8 (Asymptotics in Case I). Suppose that Assumption 4 holds and consider the regime n = dp+δ where δ ∈ (0, 1) is a constant. Fix a set S ⊆ S̄ of cardinality |S| ≤ Cdp+δ0 , with δ0 ∈ (−∞, δ), indexing polynomials of degree at most ℓ. Then for any δ ′ ∈ (0, δ − δ0 ) there exists a constant C ′ satisfying # " p ′ 1 ⊤ ≤ C ′ d−δ /2 log(n). E ΦS ΦS − I|S| n op Consequently, for any ϵ > 0, with probability at least 1 − √ C

′

log(n)

, we have ′

ΦS (X)⊤ ΦS (X)/n − I|S| op ≤ d−δ /2+ϵ . Next, we turn to the complementary setting in which n is small compared to |S|. In this case, the Gram matrix (i) (j) n ΦS Φ⊤ S = [⟨ϕS (x ), ϕS (x )⟩]i,j=1 is close to a multiple of the identity. Theorem 9 (Asymptotics in Case II). Suppose that Assumption 4 holds and consider the regime n = dp+δ where δ ∈ (0, 1) is a constant. Fix a set S ⊆ S̄ of cardinality |S| ≥ Cdp+δ0 , with δ0 ∈ (δ, ∞), indexing polynomials of degree at most ℓ. Then for any ϵ > 0 there exists a constant C ′ satisfying " # p+δ p δ0 −δ 1 E ΦS Φ⊤ ≤ C ′ d− 2 +2ϵ + log(n) · d− 2 +ϵ . (27) S − In |S| op Consequently, in the case p ≥ 1, δ0 = 1, and ϵ ∈ (0, δ), with probability at least 1 − √ C

′

log(n)

, we have

1−δ

− 2 +ϵ ΦS Φ⊤ . S /|S| − In op ≤ d

B.4

Approximating kernels by orthogonal polynomials

In this subsection, we restate several results that approximate inner-product kernels by quadratic forms in orthogonal polynomials. These results are established in [56]; we include them here for convenience, with minor notational changes, and omit the proofs. Throughout, we fix an inner-product kernel (10) satisfying the following regularity assumption. Note that Assumption 2 implies that both kernels K(z, x) and K ′ (z, x) satisfy Assumption 5. This will be important later, when we apply the results below with K replaced by K ′ . Assumption 5 (Regularity of the kernel). There exists a constant ε ∈ (0, 1) such that the function g in (10) is Lipschitz continuous on (−1 − ε, 1 + ε) and is analytic on (−ε, ε). Our goal is to approximate K(x, y) by a quadratic form of the form K(x, y) ≈ ϕ(x)⊤ Dϕ(y), where ϕ : Rd → Rp has orthonormal polynomials as its coordinate functions. Ideally, D should be relatively simple, for example diagonal. We begin by forming the Taylor approximation of g at the origin and the corresponding approximate kernel: m X g (k) (0) k ⟨x, y⟩ gm (t) := t , Km (x, y) := gm . k! d k=0

We now specialize to the uniform measure τd over the hypercube Hd = {−1, 1}d , as discussed in Section B.2. We consider the Fourier orthogonal basis ϕλ (x) := xλ for each λ ∈ {0, 1}d . As usual, for any set S ⊂ {0, 1}d , we define the concatenated vector ϕS (x) = {ϕλ (x)}λ∈S ∈ R|S| . 22

We write Φ(X) for the full matrix with entries [ϕλ (x(i) )]i,λ , and ΦS (X) for the corresponding columnsubmatrix indexed by S. It will be useful to isolate the degree-k basis elements: Sk = {λ ∈ {0, 1}d : |λ| = k}. We set ΦSk (X) and Φ≤k (X) to be, respectively, the submatrices of Φ indexed by degree k and degree at most k basis elements. We will suppress the symbol X from K(X, X) and Φ(X) throughout this subsection in order to simplify the notation. We stress that all probabilistic statements are with respect to the random vectors {x(i) }ni=1 sampled independently from τd . The next theorem, restated from [56], shows that the deviation Km (X, X) − K(X, X) is asymptotically equivalent to a multiple of the identity under favorable conditions. In particular, in the regime n = dp+δ with δ ∈ (0, 1/2), the right-hand side of (28) tends to zero for any approximation order m ≥ 2p. The approximation of kernels has been previously studied in the linear regime (i.e., n ≍ d) by [19] and in the quadratic regime (i.e., n ≍ d2 ) by [44]. Note that in [56], it is proved in the more general form for the parameterized kernel √ √ ⟨ w ⊙ x, w ⊙ y⟩ Kw (x, y) = g , d which encompasses the special case w = 1d recorded below. Theorem 10 (Taylor approximation of kernels on the hypercube). Consider independent, mean-zero, isotropic random vectors x(1) , . . . , x(n) in Rd and suppose that the coordinates of each vector x(i) are independent and uniformly distributed in {−1, 1}. Suppose that we are in the regime n = dp+δ , where δ ∈ (0, 1) is a constant. Fix an arbitrary constant C > 2 and approximation order m > 2p. Then the estimate r (C log(n))m+1 ∥K(X, X) − Km (X, X) − (g(1) − gm (1))In ∥op ≤ cg,m (28) dm−2p+1−2δ 4 holds with probability at least 1 − nC−2 , where cg,m < ∞ is a constant that depends only on m and the regularity constants in Assumption 5. ⊙k ⊤ We next record a deterministic lemma that expresses the matrix XX as a diagonal quadratic d form acting on the Fourier basis. We further decompose the diagonal matrix into the dominant part, corresponding to the degree-k basis elements, and the lower-order error part. Lemma 11 (Conversion from a polynomial kernel to a Fourier basis). For any k > 0, the following holds: ⊙k X 1 XX ⊤ e (k) Φ⊤ = d−k · ΦSk Φ⊤ ΦSj D (29) Sj , Sk + Sj k! d j: 0≤j<k, k−j is even

e (k) is a diagonal matrix with all entries on the order of Θd (d−(j+k)/2 ). where each D Sj At this point we use the permutation symmetry of the kernel (x, y) 7→ ⟨x, y⟩k . Indeed, for every permutation π of [d], one has k k ⟨x, y⟩ ⟨πx, πy⟩ = . d d Since the Walsh characters {xS yS }S⊂[d] form a basis, the coefficient of xS yS can therefore depend only on e (k) is in fact a scalar multiple of the identity, say |S|. It follows that each diagonal block D Sj

e (k) = de(k) I|S | , D j j Sj

(k) dej = Od (d−(j+k)/2 ).

Then equation (29) can be refined to ⊙k 1 XX ⊤ = d−k · ΦSk Φ⊤ Sk + k! d

X j: 0≤j<k, k−j is even

(k) dej · ΦSj Φ⊤ Sj .

(30)

The following lemma, also established in [56], is the main approximation statement in this setting. It shows the asymptotic equivalence of K(X, X) to the sum of a multiple of the identity and a quadratic form in Fourier basis elements of degree at most p. Lemma 12 (Fourier basis approximation of a kernel). Consider the regime n = dp+δ where δ ∈ (0, 1) is a constant. Then for any ϵ ∈ (0, δ), the following holds: ! log(n) ⊤ K − Φ≤p DΦ≤p − g(1) − gp (1) In op = Od,P , (31) 1−δ −ϵ 2 d where D is a diagonal matrix satisfying ∥DSk − g (k) (0)d−k I|Sk | ∥op = Od (d−k−1 ) for k = 0, . . . , p. Note that the estimate (31) allows us to approximate K by the simple expression Φ≤p DΦ⊤ ≤p + g(1) − gp (1) In up to an error of order d−(1−δ)/2 . We provide below a more refined approximation of K with the improved error rate. This sharper estimate will be useful in some of our later arguments. The proof appears in Appendix I.1. Lemma 13 (General refined Fourier-basis approximation of a kernel). Consider the regime n = dp+δ where δ ∈ (0, 1) is fixed, and let m ≥ p be a fixed integer. Then there exist scalars θp+1 , . . . , θ2m+1 , a scalar ρp,m , and a diagonal matrix D such that K − Φ≤p DΦ⊤ ≤p −

m X

⊤

θk off ΦSk ΦSk − ρp,m In

log(n)

= Od,P

k=p+1

m+1−p−δ −ϵ 2

! ,

(32)

for any sufficiently small ϵ ∈ (0, δ). Moreover, ∥DSk − g (k) (0)d−k I|Sk | ∥op = Od (d−k−1 ), θk =

(k)

(0)

k = 0, . . . , p,

+ Od (d−k−1 ),

k = p + 1, . . . , 2m + 1,

(33) (34) (35)

ρp,m = g(1) − gp (1) + Od (d−1 ). In fact one may take ρp,m := g(1) − g2m+1 (1) +

2m+1 X

θk |Sk |.

(36)

k=p+1

Proof of Lemma 1

We first fix the Walsh-degree notation used throughout the proof. Namely, for a multilinear polynomial X Y b Q(x) = Q(S) xi S⊆[d]

i∈S

and for an integer m ≥ 0, we write [Q]deg m :=

X S⊆[d]: |S|=m

b Q(S)

i∈S

for its homogeneous degree-m Walsh component. We next record the relevant subspace and the ambient basis in which coordinate estimates are expressed. Let U ∈ Rr×d have orthonormal rows (u[1] )⊤ , . . . , (u[r] )⊤ . Thus the corresponding subspace is U := row(U ) = span{u[1] , . . . , u[r] }. 24

We fix an arbitrary orthonormal completion u[r+1] , . . . , u[d] of Rd . These additional vectors are used only to represent matrices in an ambient orthonormal basis. The coherence parameter controls how much of the relevant subspace lies on each ambient coordinate. We define r d X X [j] 2 qi := ui , ρ := qi2 . j=1

i=1

Note that the coherence of U can be equivalently written as µ(U ) =

r X d [j] 2 max ui . r i∈[d] j=1

Then we have max qi = i∈[d]

rµ(U ) . d

Also, since qi = ∥PU ei ∥22 , we have 0 ≤ qi ≤ 1. Since the rows of U are orthonormal, we also have Pd i=1 qi = r. Consequently, the collision parameter satisfies ρ=

d X

d X r2 µ(U ) qi2 ≤ max qi qi = . i d i=1 i=1

(37)

We shall use the shorthand

r2 µ(U ) . d Thus all repeated-coordinate collision errors below are controlled by ρ ≤ ∆U . Notice that orthonormality alone gives ρ ≤ r, not necessarily ρ ≤ 1. The bound ρ ≤ 1 is used only in the small-coherence regime where ∆U ≤ 1. It is useful to make the finite-dimensional dependence on r explicit. We define r+ℓ r Nr,ℓ := #{λ ∈ N : |λ| ≤ ℓ} = . ℓ ∆U :=

We also define 2 2 ΘU := Nr,ℓ ∆U = Nr,ℓ

and we define 2 ΘU := rΘU = rNr,ℓ

r2 µ(U ) , d

r2 µ(U ) . d

2 The factor Nr,ℓ is a crude envelope for finite sums over multi-indices of degree at most ℓ. All constants denoted by Cℓ below depend only on ℓ, and are independent of r, d, U , and h. We now pass from the ambient variables to the relevant coordinates. We define the linear forms

zj (x) := ⟨u[j] , x⟩,

j ∈ [r].

Thus the target ridge polynomial can be written as fU∗ (x) = h(U x) = h(z1 (x), . . . , zr (x)). For simplicity, we write f ∗ instead of fU∗ when the dependence on U is clear. Qr λ For λ ∈ Nr , we write z λ := j=1 zj j . The associated all-distinct Walsh layer is defined by Hλ := H(z λ ) deg |λ| .

Hermite notation. We compare monomial and Hermite expansions of the latent polynomial h. For n ∈ N, let Hen denote the probabilists’ Hermite polynomial, namely 2

Hen (t) := (−1)n et /2

dn −t2 /2 e . dtn

(38)

Equivalently, this polynomial is given by ⌊n/2⌋

(−1)m t n−2m . m m! (n − 2m)! 2 m=0 X

Hen (t) = n! For a multi-index λ ∈ Nr , we set Heλ (z) :=

r Y

Heλj (zj ),

λ! :=

j=1

r Y

λj !.

j=1

We next introduce the index sets and coefficients appearing in the monomial-to-Hermite inversion. For α ∈ Nr , define n o A(α) := λ ∈ Nr : α − λ ∈ (2N)r . For m ≥ 0, write Am (α) := {λ ∈ A(α) : |λ| = m}. We also define

Λm (α) := {λ ∈ Nr : |λ| = m, 0 ≤ λj ≤ αj for all j ∈ [r]}.

To avoid conflict with the coherence notation, we denote the pair-count multi-index by ν(α, λ) :=

α−λ . 2

Thus, when λ ∈ A(α), the vector ν = ν(α, λ) belongs to Nr . In this case, define Aα,λ := and define Bα,λ :=

r Y

αj ! , λ ! 2νj j j=1

r Y

αj ! . λ ! 2νj νj ! j j=1

Equivalently, we have Bα,λ =

Aα,λ , ν!

ν! :=

r Y

νj !.

j=1

Applying the univariate Hermite inversion formula coordinatewise gives X zα = Bα,λ Heλ (z).

(39)

λ∈A(α)

Consequently, every polynomial h of degree at most ℓ has the unique Hermite expansion X h(z) = aλ Heλ (z). λ∈Nr : |λ|≤ℓ

We write its homogeneous Hermite components as X hq (z) := aλ Heλ (z), λ∈Nr : |λ|=q

0 ≤ q ≤ ℓ.

(40)

Similarly, we write its low-degree truncations as h≤L (z) :=

L X

0 ≤ L ≤ ℓ.

hq (z),

q=0

The corresponding latent Gaussian gradient covariance is ΣL := Ez∼γr ∇h≤L (z)∇h≤L (z)⊤ . It remains to relate the Hermite coefficients to the ordinary monomial coefficients. Write X h(z) = bα z α , Λb := {α ∈ Nr : |α| ≤ ℓ, bα ̸= 0}. |α|≤ℓ

Comparing this expansion with (39) gives X

aλ =

(41)

bα Bα,λ .

α∈Λb λ∈A(α)

We will also use the following Gaussian L2 identities. Since U has orthonormal rows, if G ∼ N (0, Id ), then U G ∼ N (0, Ir ). Therefore, ∥f ∗ ∥L2 (γd ) = ∥h∥L2 (γr ) . (42) By Hermite orthogonality under γr , we also have X

∥h∥2L2 (γr ) =

(43)

a2λ λ!.

λ∈Nr : |λ|≤ℓ

Finally, we introduce the coefficient notation used to count pair-supports. For ν ∈ Nr and T ⊆ [d], define r Y X [j] 2 Cν (T ) := [wν ] 1+ ui wj . j=1

i∈T

For the full coordinate set, we write Cν := Cν ([d]). Here w = (w1 , . . . , wr ) is a vector of formal variables. The surrogate and the target components. We now define the all-distinct surrogate. For q ∈ {0, . . . , ℓ}, set X Fq (x) := aλ Hλ (x). λ∈Nr : |λ|=q

For 0 ≤ L ≤ ℓ, define F≤L (x) :=

L X

Fq (x),

F (x) := F≤ℓ (x).

q=0

Thus the full surrogate is X

F (x) =

aλ Hλ (x).

λ∈Nr : |λ|≤ℓ

For the multilinearized target, write

P := H(f ∗ ).

We denote its Walsh-degree components by Pq := [P ]deg q ,

P≤L :=

L X

Pq ,

0 ≤ L ≤ ℓ.

q=0

The associated hypercube population matrix is M≤L := Ex∼τd ∇P≤L (x)∇P≤L (x)⊤ . Unless a measure is explicitly displayed, all L2 -norms in the hypercube part of the proof are taken with respect to x ∼ τd . 27

Proof outline. The proof proceeds in six steps. Step 1 approximates the multilinearized target P by the all-distinct surrogate F . Step 2 records the structural properties of the layers Hλ . Step 3 identifies the latent Gaussian covariance at each Hermite degree. Step 4 computes the hypercube gradient covariance of the surrogate Fq . Step 5 transfers this comparison from Fq to the true Walsh component Pq . Step 6 sums over degrees q ≤ L and obtains the covariance transfer, which completes the proof for Lemma 1. Step 1: reduction to all-distinct layers. We begin with a degree-by-degree reduction for monomials. The proof appears in Appendix I.2. Lemma 14 (Degree-m reduction). Let α ∈ Nr satisfy |α| ≤ ℓ, and let m ≤ |α| with |α| − m ∈ 2N. Then [H(z α )]deg m −

≤ Cℓ ∆U .

Aα,λ Cν(α,λ) Hλ

λ∈Am (α)

Summing over the finitely many admissible Walsh degrees gives the corresponding global reduction. Corollary 15 (Monomial reduction to all-distinct layers). For every α ∈ Nr with |α| ≤ ℓ, we have X

H(z α ) −

≤ Cℓ ∆U .

Aα,λ Cν(α,λ) Hλ

λ∈A(α)

Proof of Corollary 15. We decompose H(z α ) into homogeneous Walsh layers: X H(z α ) = [H(z α )]deg m . 0≤m≤|α| m≡|α| (mod 2)

Applying Lemma 14 to each admissible m gives an Oℓ (∆U ) error in each admissible Walsh degree. Since |α| ≤ ℓ, there are only Cℓ admissible degrees. Therefore the total error is bounded by Cℓ ∆U . This proves the claim. Next, the following lemma compares the combinatorial coefficients Cν with the Gaussian coefficients ν!−1 . The proof appears in Appendix I.3. Lemma 16 (Pair-support coefficients). For every ν ∈ Nr with |ν| ≤ ℓ, we have Cν −

1 ≤ Cℓ ∆U . ν!

The next theorem combines Corollary 15, Lemma 16, and the Hermite expansion of h then we identify the correct surrogate coefficients. The proof appears in Appendix I.4. Theorem 17 (Hypercube multilinearization in Hermite form). The following bound holds: ∥P − F ∥L2 (τd ) ≤ Cℓ ∥f ∗ ∥L2 (γd ) ΘU . Since P −F =

ℓ X

(Pq − Fq ),

q=0

and since the summands lie in distinct Walsh degrees, orthogonality gives ℓ X

∥Pq − Fq ∥2L2 = ∥P − F ∥2L2 .

q=0

(44)

Consequently, Theorem 17 implies ℓ X

∥Pq − Fq ∥2L2 ≤ Cℓ ∥f ∗ ∥2L2 (γd ) Θ2U .

(45)

q=0

In particular, for every 0 ≤ q ≤ ℓ, we have ∥Pq − Fq ∥L2 ≤ Cℓ ∥f ∗ ∥L2 (γd ) ΘU .

(46)

The same orthogonality also yields, for every cutoff 0 ≤ L ≤ ℓ, ∥P≤L − F≤L ∥2L2 =

L X

∥Pq − Fq ∥2L2 ≤ Cℓ ∥f ∗ ∥2L2 (γd ) Θ2U .

q=0

Step 2: structural properties of the all-distinct layers. We next record the structural facts about Hλ used in the covariance comparison. The proof for Proposition 18 and 19 appears in Appendix I.5 and I.6, respectively. Proposition 18 (Approximate orthogonality). Let α, β ∈ Nr satisfy |α| = |β| = m ≤ ℓ. Then ⟨Hα , Hβ ⟩ − 1{α=β} α! ≤ Cℓ ∆U . Proposition 19 (Ambient-direction contraction). Let λ ∈ Nr with 1 ≤ |λ| = q ≤ ℓ. Then, for every s ∈ [r], we have (u[s] )⊤ ∇Hλ = λs Hλ−es + Rs,λ , ∥Rs,λ ∥L2 ≤ Cℓ ∆U . Here λs Hλ−es is interpreted as 0 when λs = 0. Moreover, for every unit vector v ∈ U ⊥ , we have v ⊤ ∇Hλ = Rv,λ ,

∥Rv,λ ∥L2 ≤ Cℓ ∆U .

In particular, the second estimate applies to every vector u[s] , s ∈ [d] \ [r], in the chosen orthonormal completion. The following lemma gives a normalization comparison between the Gaussian and hypercube L2 -norms of the ridge polynomial. The proof appears in Appendix I.7. Lemma 20 (Gaussian and hypercube L2 -norm comparison). If ΘU ≤ 1, then ∥f ∗ ∥2L2 (τd ) − ∥f ∗ ∥2L2 (γd ) ≤ Cℓ ΘU ∥f ∗ ∥2L2 (γd ) . In particular, if ΘU is sufficiently small, then ∥f ∗ ∥2L2 (γd ) ≍ℓ ∥f ∗ ∥2L2 (τd ) . Step 3: the latent Gaussian covariance at degree q. For q ∈ {1, . . . , ℓ} and s, t ∈ [r], define (Gq )s,t := Eg∼N (0,Ir ) ∂s hq (g) ∂t hq (g) . (47) We extend Gq to a d × d matrix in the ambient basis {u[1] , . . . , u[d] } by setting ( (Gq )s,t , s, t ∈ [r], e (Gq )s,t := 0, otherwise. We also set Aq :=

X λ∈Nr : |λ|=q

a2λ .

Since λ! ≥ 1, Hermite orthogonality gives Aq ≤ ∥hq ∥2L2 (γr ) ≤ ∥h∥2L2 (γr ) = ∥f ∗ ∥2L2 (γd ) . By the chain rule we have

(u[s] )⊤ ∇x hq (U x) = ∂s hq (U x).

Since U x ∼ N (0, Ir ) whenever x ∼ N (0, Id ), we can derive (Gq )s,t = Ex∼N (0,Id ) (u[s] )⊤ ∇x hq (U x) (∇x hq (U x))⊤ u[t] . For 0 ≤ L ≤ ℓ, define

(48)

(49)

G≤L := Eg∼N (0,Ir ) ∇h≤L (g)∇h≤L (g)⊤ ,

which is exactly the latent Gaussian covariance ΣL . Since ∂s hq has Hermite degree q − 1, Gaussian orthogonality gives (G≤L )s,t =

L X (Gq )s,t ,

e ≤L )s,t = (G

q=1

L X

e q )s,t . (G

(50)

q=1

e ≤L denotes the extension of G≤L by zeros to the last d − r ambient coordinates. Here G Step 4: hypercube gradient covariance of the degree-q surrogate. We now compare the degree-q hypercube gradient covariance of Fq with the latent Gaussian covariance Gq . The proof appears in Appendix I.8. Proposition 21 (Ambient-basis covariance of the degree-q surrogate). Fix q ∈ {1, . . . , ℓ}, and assume ΘU ≤ 1. For s, t ∈ [d], define (Ξq )s,t := Ex∼τd (u[s] )⊤ ∇Fq (x) (∇Fq (x))⊤ u[t] . Then e q )s,t ≤ Cℓ Aq ΘU 1{s≤r or t≤r} + Θ2 1{s>r, t>r} . (Ξq )s,t − (G U The same bounds hold uniformly if either inactive basis vector u[s] , s > r, or u[t] , t > r, is replaced by an arbitrary unit vector in U ⊥ . Step 5: degreewise covariance of the multilinearized target. The next proposition transfers the degree-q comparison from Fq to Pq . The proof appears in Appendix I.9. Proposition 22 (Degreewise hypercube gradient covariance of the target). For q ∈ {1, . . . , ℓ}, define (Γq )s,t := Ex∼τd (u[s] )⊤ ∇Pq (x) (∇Pq (x))⊤ u[t] , s, t ∈ [d]. If ΘU is sufficiently small, then 2 e q )s,t ≤ Cℓ ∥f ∗ ∥2 (Γq )s,t − (G L2 (τd ) ΘU 1{s≤r or t≤r} + ΘU 1{s>r, t>r} . The same bounds hold uniformly if inactive basis vectors are replaced by arbitrary unit vectors in U ⊥ .

Step 6: covariance of the low-degree multilinearized target. We now sum over degrees. This gives the desired comparison between the hypercube gradient covariance of P≤L and the lifted Gaussian covariance of h≤L . The proof appears in Appendix I.10. Theorem 23 (Hypercube gradient covariance of the multilinearized target). For 0 ≤ L ≤ ℓ, define (Γ≤L )s,t := Ex∼τd (u[s] )⊤ ∇P≤L (x) (∇P≤L (x))⊤ u[t] , s, t ∈ [d]. If ΘU is sufficiently small, then 2 e ≤L )s,t ≤ Cℓ ∥f ∗ ∥2 (Γ≤L )s,t − (G L2 (τd ) ΘU 1{s≤r or t≤r} + ΘU 1{s>r, t>r} . The same mixed and inactive-inactive bounds hold uniformly if inactive basis vectors are replaced by arbitrary unit vectors in U ⊥ . Equivalently, since G≤L = ΣL , we have (Γ≤L )s,t − (ΣL )s,t ≤ Cℓ ∥f ∗ ∥2L2 (τd ) ΘU ,

s, t ∈ [r].

If exactly one of s, t lies in [r], then |(Γ≤L )s,t | ≤ Cℓ ∥f ∗ ∥2L2 (τd ) ΘU . Finally, for s, t ∈ [d] \ [r], we have |(Γ≤L )s,t | ≤ Cℓ ∥f ∗ ∥2L2 (τd ) Θ2U . Moreover, the corresponding block-operator comparison holds: M≤L − U ⊤ ΣL U op ≤ Cℓ ∥f ∗ ∥2L2 (τd ) ΘU .

(51)

This proves Lemma 1.

A General AGOP Approximation Result and Proof of Theorem 2

We first prove an AGOP approximation theorem that does not depend on the subspace coherence. The result is stated in coordinates aligned with the central subspace. We then apply this theorem with S = {1, . . . , r} in the rotated basis and use coherence to obtain the operator-norm bound in Theorem 2. We introduce some notation before the theorem statement. Let U⊥ ∈ R(d−r)×d have orthonormal rows spanning row(U )⊥ , and define the orthogonal completion U QU := ∈ Rd×d . U⊥ Define d×d Γq = QU Mq Q⊤ , U ∈R

Γ≤p =

Γq .

(52)

0≤q≤p

Theorem 24 (AGOP approximation). Suppose Assumptions 1, 2, and 3 hold. Suppose further that min g (k) (0) > 0,

0≤k≤p

λ+

∞ X k=p+1

g (k) (0) > 0.

Fix p ∈ {0, . . . , ℓ}, and let n = dp+δ with δ ∈ (0, 1). For every fixed deterministic subset S ⊂ [d], the following holds for any sufficiently small ϵ > 0: c − M≤p ∥op ≤ max{R(S), R(S c )} + ∥(Γ≤p )S,S c ∥op ∥M r ∥(Γ≤p )S,S ∥op + R(S) ∥(Γ≤p )S c ,S c ∥op + R(S c ) . + Here, with As := define

q=0 (Γq )s,s , Bs :=

q=0

(53)

p (Γq )s,s and AJ = (As )s∈J ∈ R|J| , BJ = (Bs )s∈J ∈ R|J| , we

p R(J) := Od,P d−δ+ϵ + d−1+δ+ϵ ∥f ∗ ∥2L2 + σε2 + Od,P d−δ/2+ϵ ∥AJ ∥1 + |J|∥AJ ∥2 p 1/2 ∥BJ ∥1 + |J|∥BJ ∥2 . + Od,P d−(1+δ)/2+ϵ/2 ∥f ∗ ∥2L2 + σε2

(54)

Proof. We establish the following key result, and defer the proof to Appendix J. Theorem 25. Suppose Assumptions 1, 2, and 3 hold. Fix p ∈ {0, . . . , ℓ}, and let n = dp+δ with δ ∈ (0, 1). Then for any ϵ > 0, and uniformly for any s, t ∈ [d], the KRR predictor fˆ satisfies h i X (u[s] )⊤ En ∇fˆ(∇fˆ)⊤ (u[t] ) − (Γq )s,t − C02 (d2δ−2 )(Γp+1 )s,t 0≤q≤p

= Od,P d

− δ2 +ϵ

# X

((Γq )s,s + (Γq )t,t ) + d

2δ−2

((Γp+1 )s,s + (Γp+1 )t,t )

0≤q≤p



ϵ − 1+δ 2 +2

+ Od,P d

X q



(Γq )s,s +

 1/2 2 (Γq )t,t  ∥f ∗ ∥L2 + σε2

0≤q≤p

q 1/2 −4+3δ ϵ q 2 (Γp+1 )s,s + (Γp+1 )t,t ∥f ∗ ∥L2 + σε2 + Od,P d 2 + 2 2 + Od,P d−1−δ+ϵ + d−2+δ+ϵ ∥f ∗ ∥L2 + σε2

(55)

where C0 = (ρ + λ)−1 g (p+1) (0). To simplify the notation, let R := ∥f ∗ ∥2L2 + σε2 ,

θd := d2δ−2 ,

κd := C02 θd .

Since QU is orthogonal, the operator norm is unchanged by this rotation, and hence c − M≤p ∥op = ∥Γ b − Γ≤p ∥op . ∥M b at the matrix appearing in Theorem 25 by defining We center Γ e := Γ b − Γ≤p − κd Γp+1 , E

b − Γ≤p . E := Γ

With this notation, the desired error is e + κd Γp+1 . E=E We first verify that the (p + 1)-degree component has bounded trace. Since Γp+1 is the rotated degree-(p + 1) population AGOP, the Fourier–Walsh gradient-energy identity gives h i 2 2 2 tr(Γp+1 ) = Ex∼τd ∥∇H(fU∗ )p+1 (x)∥2 ≲ ∥H(fU∗ )p+1 ∥L2 ≤ ∥H(fU∗ )∥L2 ≲ ∥f ∗ ∥2L2 ≤ R. 32

For the rest of the proof, introduce the temporary notation q gs := (Γp+1 )s,s , hs := (Γp+1 )s,s . By Theorem 25, uniformly over s, t ∈ [d], we have es,t | ≤ Od,P d−δ/2+ϵ [As + At + θd (gs + gt )] |E + Od,P d−(1+δ)/2+ϵ/2 R1/2 [Bs + Bt ] + Od,P d(−4+3δ)/2+ϵ/2 R1/2 [hs + ht ] + Od,P d−1−δ+ϵ + d−2+δ+ϵ R. We now convert this entrywise control into block operator-norm control. For any nonnegative vector v ∈ Rd and any J ⊂ [d], the matrix with entries vs + vt , s, t ∈ J, satisfies p ∥(vs + vt )s,t∈J ∥op ≤ ∥(vs + vt )s,t∈J ∥F ≲ ∥vJ ∥1 + |J|∥vJ ∥2 . Also, the constant entrywise term contributes at most its size times |J|, and since |J| ≤ d, it contributes Od,P d−δ+ϵ + d−1+δ+ϵ R. Applying these two deterministic bounds to the preceding entrywise estimate yields, for every deterministic J ⊂ [d], p eJ,J ∥op ≤ Od,P d−δ/2+ϵ ∥AJ ∥1 + |J|∥AJ ∥2 ∥E p + Od,P d−(1+δ)/2+ϵ/2 R1/2 ∥BJ ∥1 + |J|∥BJ ∥2 + Tp+1 (J) + Od,P d−δ+ϵ + d−1+δ+ϵ R, where Tp+1 (J) denotes the contribution of the two terms involving g and h. We next show that Tp+1 (J) is absorbed by the last line. Since Γp+1 ⪰ 0 and tr(Γp+1 ) ≲ R, we have p p √ √ ∥gJ ∥1 + |J|∥gJ ∥2 ≲ d R, ∥hJ ∥1 + |J|∥hJ ∥2 ≲ d R1/2 . Therefore the (p + 1)-degree terms satisfy Tp+1 (J) ≤ Od,P d−3(1−δ)/2+ϵ R + Od,P d−3(1−δ)/2+ϵ/2 R. Since δ ∈ (0, 1), the previous display is bounded by Tp+1 (J) ≤ Od,P d−1+δ+ϵ R. The deterministic bias κd Γp+1 is absorbed similarly. Indeed, using ∥Γp+1 ∥op ≤ tr(Γp+1 ), we get κd ∥(Γp+1 )J,J ∥op ≤ C02 d2δ−2 tr(Γp+1 ) ≲ d2δ−2 R ≤ d−1+δ R. Combining the previous bounds gives the diagonal-block estimate b − Γ≤p )J,J ∥op ≤ R(J). ∥EJ,J ∥op = ∥(Γ It remains to control the off-diagonal block. The standard PSD block inequality gives q b S,S ∥op ∥Γ b N,N ∥op . b ∥ΓS,N ∥op ≤ ∥Γ 33

Using the diagonal-block estimate with J = S and J = N , we obtain b J,J ∥op ≤ ∥(Γ≤p )J,J ∥op + R(J), ∥Γ

J ∈ {S, N }.

Consequently, the off-diagonal block of E satisfies b S,N ∥op + ∥(Γ≤p )S,N ∥op ∥ES,N ∥op ≤ ∥Γ r ∥(Γ≤p )S,S ∥op + R(S) ∥(Γ≤p )N,N ∥op + R(N ) + ∥(Γ≤p )S,N ∥op . ≤ Finally, for any symmetric block matrix one has A B ≤ max{∥A∥op , ∥D∥op } + ∥B∥op . B ⊤ D op b − Γ≤p , with the block decomposition S ∪ N = [d], gives Applying this inequality to E = Γ c − M≤p ∥op = ∥E∥op ∥M ≤ max{R(S), R(N )} + ∥(Γ≤p )S,N ∥op r + ∥(Γ≤p )S,S ∥op + R(S) ∥(Γ≤p )N,N ∥op + R(N ) . This proves the theorem.

D.1

Proof of Theorem 2

We work in an orthonormal coordinate system aligned with the relevant subspace. In this basis, N = [d] \ [r].

S = [r], Throughout the proof, r and ℓ are fixed. We write

ef := ∥f ∗ ∥2L2 + σε2 ,

µ := µ(U ).

Define the coherence scale 2 ΘU := Nr,ℓ

r2 µ , d

Nr,ℓ :=

r+ℓ . ℓ

Since r and ℓ are fixed, this scale satisfies ΘU ≲ µd . Thus, in the regime where ΘU is sufficiently small, Theorem 23 gives the block estimates ∥(Γ≤p )S,S ∥op ≲ ef , µ ∥(Γ≤p )S,N ∥op ≲ ef , d µ2 ∥(Γ≤p )N,N ∥op ≲ ef 2 . d The same theorem also gives the corresponding diagonal estimates:  ef , s ∈ S, (Γ≤p )s,s ≲ e µ2 /d2 , s ∈ N. f We next estimate the diagonal profiles that enter R(S) and R(N ). Since As takes the form As =

p X

(Γq )s,s = (Γ≤p )s,s ,

q=0

the preceding diagonal bounds imply As ≲

 ef ,

s ∈ S,

e µ2 /d2 ,

s ∈ N.

The estimate for Bs follows from Cauchy’s inequality. Indeed, because each Γq ⪰ 0, p q p X X p Bs = (Γq )s,s ≤ p + 1 (Γq )s,s q=0

!1/2 .

q=0

Since p is fixed in the present regime, this gives  e1/2 s ∈ S, f , Bs ≲ e1/2 µ/d, s ∈ N. f We now sum these coordinatewise bounds. Since |S| = r = O(1) and |N | ≤ d, we obtain p ∥AS ∥1 + |S|∥AS ∥2 ≲ ef , p 1/2 ∥BS ∥1 + |S|∥BS ∥2 ≲ ef , p µ2 ∥AN ∥1 + |N |∥AN ∥2 ≲ ef , d p 1/2 ∥BN ∥1 + |N |∥BN ∥2 ≲ ef µ. Substituting these profile estimates into the definition of R(S) gives the active-block bound R(S) ≲ ef d−δ/2+ϵ + d−1+δ+ϵ . Here the other active-block contributions are dominated by the displayed terms because δ ∈ (0, 1). Similarly, the inactive-block contribution satisfies R(N ) ≲ ef d−δ+ϵ + d−1+δ+ϵ + µ2 d−1−δ/2+ϵ + µd−(1+δ)/2+ϵ/2 . c − M≤p ∥op , the theorem yields We now apply Theorem 24. With εagop := ∥M µ εagop ≤ max{R(S), R(N )} + Cef d h i1/2 + Cef + R(S) Cef µ2 /d2 + R(N ) , where C = Cr,ℓ is a constant depending only on r and ℓ. Since ϵ > 0 is chosen sufficiently small, we have R(S) = od,P (ef ). Consequently, the active factor inside the square root is Od,P,r,ℓ (ef ). For the inactive factor, the preceding bound on R(N ) gives Cef

µ2 + R(N ) ≲ ef µ2 d−2 + d−δ+ϵ + d−1+δ+ϵ d2 + µ2 d−1−δ/2+ϵ + µd−(1+δ)/2+ϵ/2 . 35

Taking square roots, and renaming the arbitrarily small polynomial slack ϵ, the square-root term is therefore bounded by ef d−δ/2+ϵ + d−(1−δ)/2+ϵ + µd−1 + µd−1/2−δ/4+ϵ + µ1/2 d−(1+δ)/4+ϵ . Combining the preceding estimates gives εagop ≲ ef d−δ/2+ϵ + d−(1−δ)/2+ϵ + µd−1 + µ1/2 d−(1+δ)/4+ϵ + µd−1/2−δ/4+ϵ + µ2 d−1−δ/2+ϵ + µd−(1+δ)/2+ϵ/2 . Since µ ≥ 1 and δ ∈ (0, 1), the terms µd−1 and µd−(1+δ)/2+ϵ/2 are both absorbed by µd−1/2−δ/4+ϵ . Hence we conclude εagop ≲ ef d−δ/2+ϵ + d−(1−δ)/2+ϵ + µ1/2 d−(1+δ)/4+ϵ + µd−1/2−δ/4+ϵ + µ2 d−1−δ/2+ϵ . which proves (12). Finally, suppose that for some constant γ > 0, the following holds µ(U ) ≤ d1/2+δ/4−γ . Then we reach

εagop ≲ d−δ/2+ϵ + d−(1−δ)/2+ϵ + d−δ/8−γ/2+ϵ + d−γ+ϵ + d−2γ+ϵ .

Choosing ϵ > 0 sufficiently small, the right-hand side tends to zero, which completes the proof.

Proof of Lemma 3

Let ΠU := U ⊤ U and A := ΠU M≤p ΠU . We first identify the relevant eigenspace of A. Since U U ⊤ = Ir , the matrix A has the representation A = U ⊤ (U M≤p U ⊤ )U. Hence the nonzero eigenvalues of A are exactly the eigenvalues of U M≤p U ⊤ . In particular, the eigenvalue identities λr (A) = sp , λr+1 (A) = 0 hold. Moreover, the range of A is contained in row(U ), so the top-r eigenspace of A is precisely row(U ). c − A. The triangle inequality, together with the definitions of εagop Define the perturbation E := M and ρp , gives the perturbation bound c − M≤p ∥E∥op ≤ M

+ ∥M≤p − ΠU M≤p ΠU ∥op = εagop + ρp .

(56)

For brevity, set ∆ := εagop + ρp . b, U) We consider two cases. If ∆ ≥ sp /4, then 4∆/sp ≥ 1. Since sinΘ (U

≤ 1 always, the desired

estimate is immediate in this case. It remains to consider the case ∆ < sp /4. Weyl’s inequality gives the lower bound c) ≥ λr (A) − ∥E∥ ≥ sp − ∆ ≥ λ r (M op

3 sp . 4

The same inequality gives the upper bound c) ≤ λr+1 (A) + ∥E∥ ≤ ∆ ≤ λr+1 (M op

1 sp . 4

c are separated from the remaining spectrum by a These two bounds show that the top-r eigenvalues of M gap of at least sp /2. c = A + E. The r-dimensional We may now apply the Davis–Kahan sin Θ theorem to the pair A and M invariant subspace of A corresponding to its nonzero eigenvalues is row(U ), while the top-r eigenspace of c is row(U b ) by assumption. Davis–Kahan gives the estimate M b, U) sinΘ (U

≤

∥E∥op sp /2

≤2

∆ ∆ ≤4 . sp sp

Combining the two cases yields b, U) sinΘ (U

ρp + εagop ≤ min 1, 4 . sp op

Finally, under the condition (ρp + εagop )/sp = od (1), the right-hand side of (15) tends to zero. Therefore the subspace error satisfies b, U) sinΘ (U = od (1), op

b consistently recovers row(U ). which proves that U

F Set

Proof of Corollary 4 eM := M≤p − U ⊤ Σp U op ,

c − M≤p εagop := M

. op

Lemma 1 gives the deterministic comparison bound µ(U ) ∗ 2 eM = Od ∥f ∥L2 (τd ) . d The gap-stability assumption therefore implies eM = od (κ). Moreover, the weak coherence condition gives µ(U )/d = od (1), so the same bound also implies eM = od (ef ). Let PU := Prow(U ) = U ⊤ U . Since U ⊤ Σp U is supported on row(U ), the identity PU (U ⊤ Σp U )PU = U ⊤ Σp U holds. Hence the definition of ρp gives the inequality ρp = ∥M≤p − PU M≤p PU ∥op ≤ M≤p − U ⊤ Σp U op + PU (M≤p − U ⊤ Σp U )PU op ≤ 2eM . Thus ρp = od (ef ). Similarly, the compression error satisfies U M≤p U ⊤ − Σp op = U (M≤p − U ⊤ Σp U )U ⊤ op ≤ eM , where the equality uses U U ⊤ = Ir . Weyl’s inequality then yields sp = λmin (U M≤p U ⊤ ) ≥ λmin (Σp ) − eM ≥ κ − eM . Since eM = od (κ), this lower bound implies sp ≥ κ/2 for all sufficiently large d. 37

It remains to control the empirical AGOP error. Under µ(U ) = Od (d1/2+δ/4−γ ), choose ϵ > 0 sufficiently small relative to δ and γ. Then the rate Rd (U ) in Theorem 2 satisfies Rd (U ) = od (1), and therefore Theorem 2 gives the stochastic bound εagop = od,P (ef ). Combining this bound with ρp = od (ef ) gives ρp + εagop = od,P (ef ). Finally, Lemma 3 and the lower bound sp ≥ κ/2 imply b, U) sinΘ (U

≤4 op

ρp + εagop ≤ 8κ−1 (ρp + εagop ) = od,P κ−1 ef . sp

This proves the claim.

Proof of Theorem 5

In this section, we follow the notation of [37]. In particular, the space denoted by Lp elsewhere in the paper is written as Lp in this proof. Set gd := H(fU∗ ) ∈ L2 (τd ). For each k ∈ {0, . . . , d}, let Y Vd,k := span{χS : |S| = k}, χS (x) := xi . i∈S

The Walsh characters form an orthonormal basis of L2 (τd ), giving the orthogonal decomposition 2

L (τd ) =

d M

d dim(Vd,k ) = . k

Vd,k ,

k=0

Because the first-step kernel is an inner-product kernel on the hypercube, each space Vd,k is an eigenspace of 2 the associated kernel operator. Let ξd,k denote the corresponding eigenvalue, and let Pk be the orthogonal projector onto Vd,k . With this notation, gd =

d X k=0

Pk gd ,

gd,≤p :=

p X

Pk gd ,

gd,>p :=

k=0

Pk gd .

k≥p+1

By definition of the low-degree truncation of the multilinearized target, gd,≤p = H(fU∗ )≤p . We next identify the leading eigenspace appearing in Theorem 4 of [37]. Let m :=

p X d k=0

Under the assumptions of Theorem 2, the hypercube verification in Appendix D.2 of [37] applies. In particular, the degree ordering of the kernel eigenspaces is compatible with the Walsh decomposition above: the leading m eigenspaces are exactly p M

Vd,k .

k=0

Equivalently, the projector P≤m in the notation of [37] coincides with the Walsh projector P≤p := Hence P>m gd = gd,>p . 38

k=0 Pk .

2 The same hypercube spectral estimates give the fixed-degree scaling ξd,k ≍ d−k . Since p and ℓ are fixed, there exist constants c, C > 0, independent of d, such that 2 min ξd,k ≥ c d−p ,

2 max ξd,k ≤ C d−p−1 ,

0≤k≤p

p+1≤k≤ℓ

(57)

where the second bound is void when ℓ ≤ p. Let λ1 denote the regularization parameter used in the first KRR step, so that fˆ1 = fˆλ1 . By the assumptions of Theorem 2, the sample size and regularization satisfy n = dp+δ ,

λ1 ∈ [0, λ⋆ ],

0 < δ < 1,

λ⋆ := Tr(Hd,>m ).

The hypercube kernel assumptions imply Tr(Hd,>m ) = Θ(1). Therefore the effective regularization parameter from Theorem 4 of [37], γeff := λ1 + Tr(Hd,>m ), satisfies

γeff = Θ(d−p−δ ). n Theorem 4 of [37] gives the effective population estimator in diagonal form: γeff = Θ(1),

fˆγeff = eff

d X

sd,k Pk gd ,

sd,k :=

k=0

(58)

2 ξd,k . 2 ξd,k + γeff /n

For the low-degree modes k ≤ p, the lower bound in (57) and the scale estimate (58) give the shrinkage bound γeff /n γeff /n 1 − sd,k = 2 ≤ 2 ≲ d−δ . ξd,k + γeff /n ξd,k For the high-degree modes that can appear in gd , namely p + 1 ≤ k ≤ ℓ, the upper bound in (57) gives sd,k =

2 2 ξd,k ξd,k ≤ ≲ dδ−1 . 2 + γ /n ξd,k γeff /n eff

We now bound the deterministic bias of the effective estimator. Since gd has degree at most ℓ, the projection Pk gd vanishes for every k > ℓ. Orthogonality of the Walsh degree decomposition and the preceding shrinkage bounds give 2 fˆγeff − gd,≤p L2 (τ ) = eff

p X

k=0

≲ d−2δ

s2d,k ∥Pk gd ∥2L2 (τd )

k=p+1 p X

∥Pk gd ∥2L2 (τd ) + d2δ−2

k=0

≤ d

ℓ X

(1 − sd,k )2 ∥Pk gd ∥2L2 (τd ) +

−2δ

ℓ X

∥Pk gd ∥2L2 (τd )

k=p+1 2δ−2

∥gd ∥2L2 (τd ) .

(59)

Since 0 < δ < 1, both d−2δ and d2δ−2 converge to zero. Thus (59) implies 2 fˆγeff − gd,≤p L2 (τ ) = od (1) ∥gd ∥2L2 (τd ) . eff d

(60)

It remains to compare the empirical KRR estimator with the effective population estimator. Theorem 4 of [37], together with the identification P>m gd = gd,>p , gives the stochastic comparison 2 2 2 2 fˆ1 − fˆγeff = o (1) ∥g ∥ + ∥g ∥ + σ (61) 2 2+η d,P d L (τd ) d,>p L ε . (τd ) eff L2 (τ ) d

We control the L2+η -term using the bounded degree of gd . Since deg(gd ) ≤ ℓ, the function gd,>p also has degree at most ℓ. Set q0 := ⌈2 + η⌉. Monotonicity of Lq -norms on a probability space and the Boolean hypercontractivity estimate from Lemma 18 of [37] give ∥gd,>p ∥2L2+η (τd ) ≤ ∥gd,>p ∥2Lq0 (τd ) ≤ (q0 − 1)ℓ ∥gd,>p ∥2L2 (τd ) . Since gd,>p is an orthogonal projection of gd , the preceding estimate implies ∥gd,>p ∥2L2+η (τd ) ≲ ∥gd ∥2L2 (τd ) . Substituting this bound into (61) yields 2 2 2 fˆ1 − fˆγeff = o (1) ∥g ∥ + σ 2 d,P d L (τd ) ε . eff L2 (τ )

(62)

Finally, we combine the empirical-to-effective error with the deterministic effective bias. The elementary inequality ∥a + b∥22 ≤ 2∥a∥22 + 2∥b∥22 , applied in L2 (τd ) with a = fˆ1 − fˆγeff , eff

b = fˆγeff − gd,≤p , eff

gives 2

∥fˆ1 − gd,≤p ∥2L2 (τd ) ≤ 2 fˆ1 − fˆγeff + 2 fˆγeff − gd,≤p L2 (τ ) eff L2 (τd ) eff d 2 2 2 = od,P (1) ∥gd ∥L2 (τd ) + σε + od (1) ∥gd ∥L2 (τd ) = od,P (1) ∥gd ∥2L2 (τd ) + σε2 .

(63)

Since gd,≤p = H(fU∗ )≤p , the estimate (63) is the desired claim.

Proof of Proposition 6

Let

B := U ⊤ Σp U,

c1 − B ∆d := M

By construction, M2 =

md := max{d−δ/2 , d−1+δ }.

, op

d c (M1 + ηId ), cη

c1 + ηId ). cη = tr(M

We first identify the deterministic approximation term. Let U Q := ∈ O(d). U⊥ Since U⊥ spans row(U⊥ ), we have

Id = U ⊤ U + (U⊥ )⊤ U⊥ .

Using this decomposition, we may write Σp + ηIr B + ηId = U (Σp + ηIr )U + η(U⊥ ) U⊥ = Q 0 ⊤

⊤

Taking principal square roots yields p B + ηId = Q⊤

Σp + ηIr 0

√

0 Q. η Id−r

0 ηId−r

Therefore, by the definition of x b, we have x b=

B + ηId x.

This identity allows us to rewrite the approximation error as s q q p d 1/2 d c1 + ηId − B + ηId x b L = M M2 x − cη x L2 2 cη s q p d c1 + ηId − B + ηId ≤ ∥x∥L2 . M cη op c1 + ηId and B + ηId are bounded We next control the square-root perturbation term. Since both M below by ηId , the square-root map is operator-Lipschitz on [η, ∞), which gives q p ∆d 1 c1 − B c1 + ηId − B + ηId = √ . ≤ √ M M 2 η 2 η op op Substituting this bound into the previous display, we obtain √ q d 1/2 d b L2 ≤ √ ∆d ∥x∥L2 . M2 x − cη x 2 cη η It remains to identify the scale of cη . Expanding the trace, we have c1 − B). cη = ηd + tr(B) + tr(M Since U U ⊤ = Ir , this gives

tr(B) = tr(U ⊤ Σp U ) = tr(Σp ),

and therefore we have |cη − ηd| ≤ tr(Σp ) + d ∆d . Now Corollary 4, applied with ϵ = ζ/2, yields ∆d = Od,P dζ/2 md . Because η = dζ md , we have

d ∆d = Od,P d−ζ/2 ηd = od,P (ηd).

Moreover, the following holds

ηd ≥ dζ · d−1+δ · d = dζ+δ → ∞,

while tr(Σp ) = Od (1). Hence tr(Σp ) = o(ηd), and we conclude that cη = (1 + od,P (1))ηd. Finally, inserting this asymptotic into (64) yields √ d 1 = (1 + od,P (1)) , √ cη η η and therefore we obtain 1/2

M2 x −

∆d d b L2 ≤ (1 + od,P (1)) ∥x∥L2 . cη x 2η

Using again ∆d = Od,P (dζ/2 md ) and η = dζ md , we obtain ∆d = Od,P (d−ζ/2 ), η which gives 1/2

M2 x −

d b L2 = Od,P (d−ζ/2 ) ∥x∥L2 . cη x

This is exactly (20). 41

(64)

Proof for preliminaries

I.1

Proof of Lemma 13

We begin with the polynomial truncation. Applying Theorem 10 with truncation order 2m + 1, we obtain m+1 log (n) . (65) ∥K − K2m+1 − (g(1) − g2m+1 (1)) In ∥op = Od,P dm+1−p−δ We next expand K2m+1 in the Fourier basis. By definition of the truncated kernel, we have K2m+1 =

2m+1 X

g (ℓ) (0) ·

ℓ=0

1 ℓ!

XX ⊤ d

⊙ℓ .

Applying the refined identity (30) term by term, this gives  K2m+1 =

2m+1 X ℓ=0

 −ℓ ⊤ g (ℓ) (0)  d ΦSℓ ΦSℓ +

 X j: 0≤j<ℓ, ℓ−j is even

 (ℓ)  dej ΦSj Φ⊤ Sj  .

Collecting together the terms with the same Fourier degree, we may therefore write K2m+1 = Φ≤p DΦ⊤ ≤p +

2m+1 X

θk ΦSk Φ⊤ Sk ,

(66)

k=p+1

where, for each j = 0, . . . , p, 

  g (j) (0) DSj :=   dj +

X ℓ: j<ℓ≤2m+1, ℓ−j is even

(ℓ)  g (ℓ) (0) dej   I|Sj | ,

and, for each k = p + 1, . . . , 2m + 1, θk :=

g (k) (0) + dk

(ℓ)

g (ℓ) (0) dek .

ℓ: k<ℓ≤2m+1, ℓ−k is even

(ℓ) These formulas immediately yield the required coefficient bounds. Indeed, since dej = Od (d−(j+ℓ)/2 ) and ℓ − j is a positive even integer in the correction terms, we necessarily have ℓ ≥ j + 2. It follows that (ℓ) dej = Od (d−j−1 ).

Because the number of indices ℓ is bounded in terms of m, the whole correction sum is still of order Od (d−j−1 ). Hence, for j = 0, . . . , p, we obtain ∥DSj − g (j) (0)d−j I|Sj | ∥op = Od (d−j−1 ), while for k = p + 1, . . . , 2m + 1, we have θk =

g (k) (0) + Od (d−k−1 ). dk

This proves the coefficient estimates in the statement. 42

I.2

Proof of Lemma 14

Let X ∼ τd . For each S ⊆ [d] with |S| = m, set aS := E[z α (X)X S ]. The degree-m Walsh component of H(z α ) is X [H(z α )]deg m (x) = aS xS . |S|=m

For λ ∈ Λm (α), write X

Hλ (x) =

hλ,S xS ,

hλ,S := E[z λ (X)X S ].

|S|=m

Qr α Fix S ⊆ [d] with |S| = m. Expanding z α = j=1 zj j amounts to assigning each slot corresponding to direction j to a coordinate i ∈ [d]. Such an assignment contributes to aS precisely when, for each coordinate i, the total multiplicity assigned to i has parity 1{i∈S} . We first count the principal assignments. Fix λ ∈ Am (α), and put ν = ν(α, λ) = (α − λ)/2. A contributing assignment is called principal of type λ if it has the following structure: each coordinate i ∈ S is hit exactly once; among these singleton hits, exactly λj come from direction j; the remaining 2νj slots of direction j are grouped into νj same-direction pairs; and the pair-support coordinates are distinct and lie in S c . Let Gλ,S denote the total contribution of the principal assignments of type λ. A direct count gives the identity Gλ,S = Aα,λ hλ,S Cν(α,λ) (S c ). (67) The factor hλ,S accounts for the singleton placements on S, while Cν(α,λ) (S c ) encodes the pair-support coordinates in S c . For direction j, the slot combinatorics contribute the factor αj ! αj (2νj )! νj ! = . ν j λj ! 2νj λj 2 νj ! This is the j-th factor in Aα,λ =

r Y

αj ! . λ ! 2νj j=1 j

Summing over all admissible λ, the coefficient aS decomposes as X aS = Gλ,S + BS , λ∈Am (α)

where BS denotes the total contribution of the non-principal assignments. We next replace Cν (S c ) by Cν . Since all coefficients in the generating function defining Cν (T ) are nonnegative, the difference Cν − Cν (S c ) is obtained by forcing at least one pair-support coordinate to lie in S. Since |ν| ≤ ℓ, this difference satisfies X |Cν − Cν (S c )| ≤ Cℓ qi . i∈S

Because |S| = m ≤ ℓ, the preceding bound is at most Cℓ m maxi qi ≤ Cℓ ∆U . Using (67), the finiteness of Am (α), and the bound ∥Hλ ∥L2 ≤ ∥z λ ∥L2 ≤ Cℓ , the omission error satisfies 2

Aα,λ Cν(α,λ) (S c ) − Cν(α,λ) hλ,S

≤ Cℓ ∆2U .

(68)

|S|=m λ∈Am (α)

It remains to bound the non-principal assignments. We first record an elementary coefficient estimate. For a fixed ordered label list σ = (σ1 , . . . , σm ) ∈ [r]m , define Kσ,S :=

m Y

ψ: [m]→S bijection a=1

[σ ]

a uψ(a) .

The following bound will be used repeatedly: X

(69)

2 Kσ,S ≤ Cℓ .

|S|=m

Indeed, Kσ,S is the absolute-value analogue of a degree-m Walsh coefficient of a product of m linear forms. Since the degree-m projection is an L2 (τd )-contraction, Hölder’s inequality and the Khintchine inequality give  1/2 ! m d m d X Y X Y X [σa ] [σ ] 2   Kσ,S ≤ |ui |Xi ≤ |ui a |Xi ≤ Cℓ . |S|=m

a=1

i=1

a=1

i=1

L2m

This proves (69). We now decompose the non-principal assignments. We overcount them by making two boundedcomplexity choices. First, for each coordinate i ∈ S, choose one distinguished hit among the odd number of hits at i. Second, choose one distinguished witness of non-principality. Since |α| ≤ ℓ, the number of such choices is bounded by Cℓ . The distinguished singleton hits on S have an ordered label list σ = (σ1 , . . . , σm ) ∈ [r]m , and their total absolute contribution is bounded by Kσ,S . After these distinguished singleton hits are fixed, all remaining multiplicities are even. We impose the following priority rule for the remaining defects. 1. First, consider assignments for which at least one remaining block touches S. Choose such a block as the distinguished witness. Let j1 , . . . , jt , with t ≥ 2, be its fixed labels. For a support coordinate i ∈ S, this block contributes at most t Y

[j ]

t/2

|ui a | ≤ qi

≤ qi ,

a=1

where the last inequality uses qi ≤ 1. Summing over the possible support coordinate in S gives the bound X qi ≤ m max qi ≤ Cℓ ∆U . i

i∈S

2. Second, consider the remaining assignments for which no remaining block touches S, but some coordinate in S c carries remaining multiplicity at least 4. Choose such a block as the distinguished witness. If its fixed labels are j1 , . . . , jt , with t ≥ 4, then its total absolute contribution is at most d Y t X

[j ] |ui a | ≤

i=1 a=1

d d X X t/2 qi ≤ qi2 = ρ ≤ ∆U . i=1

i=1

3. Third, consider the remaining assignments for which no remaining block touches S, and no coordinate in S c carries remaining multiplicity at least 4. At this stage, every remaining block is a 2-block in S c , and the supports of these 2-blocks are distinct. If the assignment is still non-principal, then at least one such 2-block is mixed-direction. Choose a mixed block with labels j = ̸ k as the distinguished witness. After all other remaining blocks have been fixed, the allowed support coordinates for this mixed block are i ∈ / T , where T contains S and the supports of the other remaining blocks. In particular, |T | ≤ Cℓ . The orthogonality of the rows of U gives the identity X [j] [k] X [j] [k] ui ui . ui ui = − i∈T

i∈T /

This identity bounds the absolute contribution of the mixed block by X i∈T /

[j] [k]

ui ui

≤

qi ≤ Cℓ max qi ≤ Cℓ ∆U .

i∈T

It remains to control the non-distinguished even blocks. For a fixed block with labels j1 , . . . , jt , where t ≥ 2, Hölder’s inequality gives d Y t X

[j ]

|ui a | ≤

i=1 a=1

t Y

∥u[ja ] ∥t ≤

a=1

∥u[ja ] ∥2 = 1.

a=1

Restrictions on the allowed support coordinates can only decrease this absolute sum. Since the number of remaining blocks is bounded by ℓ, all non-distinguished even blocks together contribute at most Cℓ . The three prioritized defect estimates give the pointwise bound X |BS | ≤ Cℓ ∆U Kσ,S , σ∈Sα,m

where Sα,m is a finite set of ordered label lists with cardinality bounded by Cℓ . Squaring, summing over S, and using (69) gives X |BS |2 ≤ Cℓ ∆2U . (70) |S|=m

Finally, define the coefficient error by X

ES := aS −

Aα,λ Cν(α,λ) hλ,S .

λ∈Am (α)

The estimates (68) and (70) imply X

|ES |2 ≤ Cℓ ∆2U .

|S|=m

By Parseval’s identity, the left-hand side is exactly the squared L2 -error of the degree-m approximation. Taking square roots proves the lemma.

I.3

Proof of Lemma 16

Define Li (w) := d Y

[j] 2 j=1 (ui ) wj . The generating function factorizes as

(1 + Li (w)) = exp

i=1

d X

! Li (w) exp(R(w)),

i=1 [j]

[wν ] exp

log(1 + Li (w)) − Li (w) .

i=1

Since i (ui )2 = 1 for every j, the linear term satisfies wν in the first exponential factor is P

R(w) :=

d X

! Li (w)

i=1

i=1 Li (w) =

j=1 wj . Thus the coefficient of

  r X 1 = [wν ] exp  wj  = . ν! j=1

We first consider the case ∆U ≥ 1. Since 1 + Li (w) ≤ exp(Li (w)) coefficientwise, the coefficient Cν satisfies ! d X 1 ν 0 ≤ Cν ≤ [w ] exp Li (w) = . ν! i=1 In this case, the desired error bound follows from Cν −

1 1 ≤ ≤ Cℓ ∆U . ν! ν!

It remains to consider the case ∆U < 1. In this case, ρ ≤ ∆U < 1. The series R(w) has no constant or linear terms. Therefore, every coefficient of R(w) up to total degree |ν| ≤ ℓ is a finite linear combination, 45

with coefficients depending only on ℓ, of coefficients of Li (w)m , where 2 ≤ m ≤ |ν|. For such m, the coefficient bound |[wω ]Li (w)m | ≤ Cℓ qim ≤ Cℓ qi2 holds because 0 ≤ qi ≤ 1. After summing over i, every relevant coefficient of R(w) is bounded by Cℓ

d X

qi2 = Cℓ ρ.

i=1

Since ρ< 1 and|ν| ≤ ℓ, every coefficient of exp(R(w))−1 up to total degree |ν| is also Oℓ (ρ). Multiplication P by exp j wj , whose coefficients up to degree ℓ are bounded by Cℓ , gives Cν −

1 ≤ Cℓ ρ ≤ Cℓ ∆U . ν!

This proves the lemma.

I.4

Proof of Theorem 17

Write the monomial expansion of h as h(z) =

bα z α .

α∈Λb

We first record the coefficient comparison X X 1/2 |bα | + |aλ | ≤ Cℓ Nr,ℓ ∥h∥L2 (γr ) . |α|≤ℓ

(71)

|λ|≤ℓ

Indeed, Hermite orthogonality gives 1/2

 X

1/2

|aλ | ≤ Nr,ℓ 

|λ|≤ℓ

Using (43), the preceding display becomes X

a2λ λ!

|λ|≤ℓ

1/2

|aλ | ≤ Nr,ℓ ∥h∥L2 (γr ) .

|λ|≤ℓ

The same bound for the monomial coefficients follows by expanding each Heλ into monomials. Since |λ| ≤ ℓ, the ℓ1 -norm of the monomial coefficient vector of Heλ is bounded by Cℓ . This proves (71). Applying Corollary 15 to each monomial in h, and then using (71), gives X H(f ∗ ) = cλ Hλ + E, |λ|≤ℓ

where

1/2

∥E∥L2 (τd ) ≤ Cℓ Nr,ℓ ∥h∥L2 (γr ) ∆U . Here the coefficients cλ are given by cλ =

bα Aα,λ Cν(α,λ) .

α∈Λb λ∈A(α)

On the other hand, the Hermite coefficient formula (41) gives aλ =

bα Aα,λ

α∈Λb λ∈A(α)

1 . ν(α, λ)!

Subtracting the two coefficient formulas and applying Lemma 16 yields 1/2

|cλ − aλ | ≤ Cℓ Nr,ℓ ∥h∥L2 (γr ) ∆U . Moreover, the layer norm satisfies ∥Hλ ∥L2 (τd ) ≤ ∥z λ ∥L2 (τd ) ≤ Cℓ , where the last inequality follows from Hölder’s inequality and the Khintchine inequality. Since there are at most Nr,ℓ indices λ, we get X

3/2

(cλ − aλ )Hλ

|λ|≤ℓ

≤ Cℓ Nr,ℓ ∥h∥L2 (γr ) ∆U . L2 (τd )

3/2

2 Since Nr,ℓ ≤ Nr,ℓ , the preceding display is bounded by Cℓ ∥h∥L2 (γr ) ΘU . Combining this estimate with the bound on E, and using (42), proves

∥P − F ∥L2 (τd ) ≤ Cℓ ∥f ∗ ∥L2 (γd ) ΘU . This completes the proof.

I.5

Proof of Proposition 18

The cases m = 0 and m = 1 are immediate, so assume m ≥ 2. Choose ordered label lists σ1 , . . . , σm ∈ [r] and τ1 , . . . , τm ∈ [r] such that #{p : σp = j} = αj ,

#{p : τp = j} = βj .

With this choice, the two monomials can be written as zα =

m Y

zβ =

zσp ,

p=1

m Y

zτp .

p=1

Expanding each linear form shows that the top-degree layers keep exactly the injective terms. Thus X

Hα (x) =

m Y

[σ ] uip p xi1 · · · xim ,

i1 ,...,im ∈[d] p=1 all distinct

and similarly X

Hβ (x) =

m Y

p] u[τ x v1 · · · x vm . vp

v1 ,...,vm ∈[d] p=1 all distinct

Let X ∼ τd . The expectation is nonzero precisely when the two injective tuples have the same underlying coordinate set. Hence we have ⟨Hα , Hβ ⟩ =

m X Y

i1 ,...,im π∈Sm p=1 all distinct

[σ ] [τ −1 (p) ]

uip p uipπ

We compare this restricted sum with the unrestricted sum m X Y

M :=

[σ ] [τ −1 (p) ]

uip p uipπ

i1 ,...,im ∈[d] π∈Sm p=1

For fixed π, the sum over (i1 , . . . , im ) factorizes as m X Y

[σ ] [τ −1 (p) ]

uip p uipπ

m Y

i1 ,...,im p=1

⟨u[σp ] , u[τπ−1 (p) ] ⟩.

p=1

By orthonormality, this product equals 1 exactly when σp = τπ−1 (p) for every p, and otherwise it equals 0. Such permutations exist if and only if α = β. In that case, their number is r Y

αj ! = α!.

j=1

Therefore, we get M = 1{α=β} α!. It remains to control the contribution of tuples with collisions. Since ⟨Hα , Hβ ⟩ is obtained from M by restricting to tuples with all entries distinct, we have M − ⟨Hα , Hβ ⟩ =

X π∈Sm

m Y

[σ ] [τ −1 (p) ]

uip p uipπ

p=1 i1 ,...,im not all distinct

Using the elementary bound X

1{not all distinct} ≤

1{ip =iq } ,

1≤p<q≤m

and then taking absolute values, we get X

|M − ⟨Hα , Hβ ⟩| ≤

m X Y

[σ ] [τ −1 (c) ]

uic c uicπ

π∈Sm 1≤p<q≤m i1 ,...,im c=1 ip =iq

Fix π ∈ Sm and p < q, and write i = ip = iq . For any j, k ∈ [r], we have [j] [k]

|ui ui | ≤

1 [j] [k] (ui )2 + (ui )2 ≤ qi . 2

Thus the two constrained positions contribute at most qi2 . Pd [σ ] [τ −1 ] For each remaining position, Cauchy–Schwarz gives v=1 |uv c uv π (c) | ≤ 1. Therefore each conPd strained collision sum is bounded by i=1 qi2 = ρ. Summing over π ∈ Sm and over all pairs p < q, we obtain m ⟨Hα , Hβ ⟩ − 1{α=β} α! ≤ m! ρ. 2 Since m ≤ ℓ and ρ ≤ ∆U , the desired bound follows.

I.6

Proof of Proposition 19

We first prove the active-direction estimate. Fix s ∈ [r], and choose an ordered label list σ1 , . . . , σq ∈ [r] such that #{p : σp = j} = λj for each j ∈ [r]. The all-distinct layer Hλ admits the injective expansion ! q X Y [σp ] Hλ (x) = utp xt1 · · · xtq . t1 ,...,tq ∈[d] all distinct

p=1

For a tuple (tv )v̸=p , write S = {tv : v = ̸ p}. Differentiating the preceding expansion term by term gives the directional derivative formula   q X X Y [σ ] X [s] [σ ]  (u[s] )⊤ ∇Hλ (x) = utvv  xS ut ut p . p=1

(tv )v̸=p all distinct

t∈S /

v̸=p

The orthonormality of the rows of U gives, for each j ∈ [r], the identity X [s] [j] X [s] [j] ut ut = 1{s=j} − ut ut . t∈S

t∈S /

Substituting this identity into the derivative formula separates the main term from the error term. The contribution of 1{s=σp } is exactly λs Hλ−es , with the convention that this term is zero when λs = 0. The remaining terms define Rs,λ , giving the decomposition (u[s] )⊤ ∇Hλ = λs Hλ−es + Rs,λ . It remains to estimate the remainder. For every (q − 1)-element set S, the correction factor satisfies X

[s] [σ ]

ut ut p ≤ |S| max qi ≤ Cℓ ∆U . i

t∈S

After grouping the p-th summand by the underlying set S, its Walsh coefficients are those of Hλ−eσp multiplied by this correction factor. Since ∥Hλ−eσp ∥L2 ≤ ∥z λ−eσp ∥L2 ≤ Cℓ , Parseval’s identity gives (p)

∥Rs,λ ∥L2 ≤ Cℓ ∆U . Summing over p ∈ [q], with q ≤ ℓ, yields ∥Rs,λ ∥L2 ≤ Cℓ ∆U . We now prove the inactive-direction estimate. Let v ∈ U ⊥ be a unit vector. The same injective expansion gives   q X X Y [σ ] X [σ ] a ⊤  uta  xS vt ut p . v ∇Hλ (x) = p=1

(ta )a̸=p all distinct

a̸=p

[σ ]

t∈S /

Since v ⊥ u[σp ] , the inner sum satisfies vt ut p = −

[σ ]

vt ut p .

t∈S

t∈S /

Thus v ⊤ ∇Hλ consists only of remainder terms; write v ⊤ ∇Hλ = Rv,λ . (p) Let cS denote the coefficient of xS in Hλ−eσp . By Cauchy–Schwarz, the coefficient multiplying xS in the p-th remainder satisfies 2 X X [σp ] [σ ] vi ui ≤ (q − 1) vi2 (ui p )2 . i∈S

i∈S

Consequently, Parseval’s identity gives the preliminary bound (p)

∥Rv,λ ∥2L2 ≤ Cℓ

d X

[σ ]

vi2 (ui p )2

i=1

X S∋i

(p)

|cS |2 .

(72)

Claim 1. For every i ∈ [d], the coefficient localization bound holds: X (p) |cS |2 ≤ Cℓ qi . S∋i (p)

Proof of Claim 1. Set η := λ − eσp . Since cS is the coefficient of xS in Hη , Parseval’s identity gives X (p) |cS |2 = ∥∂i Hη ∥2L2 . S∋i

Choose an ordered label list ρ1 , . . . , ρq−1 for η. The injective expansion of Hη is ! q−1 Y [ρ ] X a uta xt1 · · · xtq−1 . Hη (x) = a=1

t1 ,...,tq−1 ∈[d] all distinct

Differentiating this expansion with respect to xi exposes one label at a time. In the term exposing the [ρ ] label ρa , the coefficient ui a appears, and the remaining factor is an all-distinct layer of degree q − 2 with coordinate i excluded from its support. The same Khintchine and Hölder bounds used above show that this remaining layer has L2 (τd )-norm at most Cℓ . Therefore the derivative norm satisfies ∥∂i Hη ∥L2 ≤ Cℓ

q−1 X

[ρ ]

|ui a |.

a=1

Since q ≤ ℓ, Cauchy–Schwarz gives the bound ∥∂i Hη ∥2L2 ≤ Cℓ

q−1 X

[ρ ]

(ui a )2 .

a=1

The exposed labels are among the relevant directions, so q−1 X

[ρ ]

(ui a )2 ≤ Cℓ

a=1

r X

[j]

(ui )2 = Cℓ qi .

j=1

Combining the preceding two estimates proves the claim. Substituting Claim 1 into (72) gives (p)

∥Rv,λ ∥2L2 ≤ Cℓ

d X

[σ ]

vi2 (ui p )2 qi .

i=1 [σ ]

Since (ui p )2 ≤ qi , the last sum satisfies d X

[σ ]

vi2 (ui p )2 qi ≤

i=1

Using maxi qi ≤ ∆U , we obtain

d X

vi2 qi2 ≤ (max qi )2 .

i=1

(p)

∥Rv,λ ∥2L2 ≤ Cℓ ∆2U . Equivalently, the p-th remainder satisfies (p)

∥Rv,λ ∥L2 ≤ Cℓ ∆U . Summing over p ∈ [q], with q ≤ ℓ, gives ∥Rv,λ ∥L2 ≤ Cℓ ∆U . This completes the proof. 50

I.7

Proof of Lemma 20

Since P = H(f ∗ ) agrees with f ∗ on {±1}d , we have ∥P ∥L2 (τd ) = ∥f ∗ ∥L2 (τd ) . By Theorem 17, there is a decomposition P = F + E, ∥E∥L2 (τd ) ≤ Cℓ ∥f ∗ ∥L2 (γd ) ΘU . (73) Using (42), the error bound in (73) may equivalently be written as ∥E∥L2 (τd ) ≤ Cℓ ∥h∥L2 (γr ) ΘU . We first compare the surrogate F with the latent polynomial h. Different Walsh degrees are exactly orthogonal, and Proposition 18 controls the inner products within each fixed degree. Therefore the squared norms satisfy  2 X X ∥F ∥2L2 (τd ) − a2λ λ! ≤ Cℓ ∆U  |aλ | . |λ|≤ℓ

|λ|≤ℓ

2 2 |λ|≤ℓ aλ λ! = ∥h∥L2 (γr ) , the coefficient bound (71), and the inequality Nr,ℓ ∆U ≤ ΘU

The Hermite identity give the key comparison P

∥F ∥2L2 (τd ) − ∥h∥2L2 (γr ) ≤ Cℓ ΘU ∥h∥2L2 (γr ) .

(74)

Since ΘU ≤ 1, the comparison (74) also implies ∥F ∥L2 (τd ) ≤ Cℓ ∥h∥L2 (γr ) . We now transfer the comparison from F to P . The decomposition (73) gives the deterministic bound ∥P ∥2L2 (τd ) − ∥F ∥2L2 (τd ) ≤ 2∥F ∥L2 (τd ) ∥E∥L2 (τd ) + ∥E∥2L2 (τd ) ≤ Cℓ ΘU ∥h∥2L2 (γr ) .

(75)

In the last step, we used the preceding bound on ∥F ∥L2 (τd ) , the error estimate for E, and the assumption ΘU ≤ 1. Combining (74) and (75) gives ∥P ∥2L2 (τd ) − ∥h∥2L2 (γr ) ≤ Cℓ ΘU ∥h∥2L2 (γr ) . Finally, using ∥P ∥L2 (τd ) = ∥f ∗ ∥L2 (τd ) and (42), we obtain the desired estimate ∥f ∗ ∥2L2 (τd ) − ∥f ∗ ∥2L2 (γd ) ≤ Cℓ ΘU ∥f ∗ ∥2L2 (γd ) . If ΘU is sufficiently small, then (76) implies the two-sided bound (1 − Cℓ ΘU )∥f ∗ ∥2L2 (γd ) ≤ ∥f ∗ ∥2L2 (τd ) ≤ (1 + Cℓ ΘU )∥f ∗ ∥2L2 (γd ) . This proves the stated norm equivalence.

I.8

Proof of Proposition 21

For s ∈ [d], define

Ds := (u[s] )⊤ ∇Fq .

By Proposition 19, we have

Ds =

X  aλ λs Hλ−es + Rs , 

s ∈ [r],

|λ|=q

 R ,

s ∈ [d] \ [r],

(76)

where the remainder satisfies

∥Rs ∥L2 ≤ Cℓ ∆U

|aλ |.

|λ|=q

Nr,q := #{λ ∈ Nr : |λ| = q},

then Cauchy–Schwarz gives X

1/2 1/2 |aλ | ≤ Nr,q Aq .

|λ|=q

Therefore, we have

1/2 ∥Rs ∥L2 ≤ Cℓ Nr,q ∆U A1/2 q .

Assume first that s ∈ [r]. Reindexing with γ = λ − es gives X Ds = (γs + 1)aγ+es Hγ + Rs . |γ|=q−1

Denote the leading term by X

Ls :=

(γs + 1)aγ+es Hγ .

|γ|=q−1

Using Proposition 18 within degree q − 1, we obtain  ∥Ls ∥2L2 ≤ Cℓ Aq + Cℓ ∆U 

2 X

|(γs + 1)aγ+es | .

|γ|=q−1

By Cauchy–Schwarz, this implies ∥Ls ∥2L2 ≤ Cℓ Aq + Cℓ Nr,q ∆U Aq . 1/2

2 Since Nr,q ≤ Nr,ℓ and ΘU = Nr,ℓ ∆U ≤ 1, we get ∥Ls ∥L2 ≤ Cℓ Aq . Together with the remainder bound, this yields ∥Ds ∥L2 ≤ Cℓ A1/2 s ∈ [r]. q ,

If s ∈ [d] \ [r], then Ds = Rs , and hence ∥Ds ∥L2 ≤ Cℓ A1/2 q ΘU . We first consider the active-active block. If s, t ∈ [r], then expanding the covariance gives E[Ds Dt ] = ⟨Ls , Lt ⟩ + Oℓ (ΘU Aq ). Using Proposition 18, the leading inner product satisfies ⟨Ls , Lt ⟩ = (Gq )s,t + Oℓ (Nr,q ∆U Aq ). Since Nr,q ∆U ≤ ΘU , we obtain e q )s,t ≤ Cℓ ΘU Aq , (Ξq )s,t − (G

s, t ∈ [r].

Next suppose exactly one of s, t is inactive. Then we immediately have e q )s,t = 0. (G 1/2

1/2

One derivative has L2 -norm at most Cℓ Aq , while the other has L2 -norm at most Cℓ Aq ΘU . Therefore, we have e q )s,t ≤ Cℓ ΘU Aq . (Ξq )s,t − (G 52

Finally, if s, t > r, then both derivatives are inactive remainders. Thus we have e q )s,t ≤ Cℓ Nr,q ∆2U Aq . (Ξq )s,t − (G Since Nr,q ≤ Nr,ℓ , this is bounded by Cℓ Θ2U Aq . Combining the three cases proves the entrywise estimate. The same proof works uniformly if any inactive basis vector is replaced by an arbitrary unit vector in U ⊥ , because Proposition 19 is uniform over such unit vectors.

I.9

Proof of Proposition 22

Set Eq := Pq − Fq . By (46), we have ∥Eq ∥L2 ≤ Cℓ ∥f ∗ ∥L2 (γd ) ΘU . Since Eq is homogeneous multilinear of degree q, for any v ∈ Rd , its directional derivative has the expansion   X X cq (T ∪ {j}) xT .  v ⊤ ∇Eq = vj E |T |=q−1

By Cauchy–Schwarz, this implies

j ∈T /

∥v ⊤ ∇Eq ∥2L2 ≤ q∥v∥22 ∥Eq ∥2L2 .

Thus, for every unit vector v, we have ∥v ⊤ ∇Eq ∥L2 ≤ Cℓ ∥f ∗ ∥L2 (γd ) ΘU . For convenience, write P Ds,q := (u[s] )⊤ ∇Pq ,

F Ds,q := (u[s] )⊤ ∇Fq ,

E Ds,q := (u[s] )⊤ ∇Eq .

Since Pq = Fq + Eq , these derivatives satisfy P F E Ds,q = Ds,q + Ds,q .

Expanding the covariance difference gives E F F E E E (Γq )s,t − (Ξq )s,t = E[Ds,q Dt,q ] + E[Ds,q Dt,q ] + E[Ds,q Dt,q ].

Taking absolute values and applying Cauchy–Schwarz yields E F F E E E (Γq )s,t − (Ξq )s,t ≤ ∥Ds,q ∥L2 ∥Dt,q ∥L2 + ∥Ds,q ∥L2 ∥Dt,q ∥L2 + ∥Ds,q ∥L2 ∥Dt,q ∥L2 .

The derivative bound for Eq gives E ∥Ds,q ∥L2 ≤ Cℓ ∥f ∗ ∥L2 (γd ) ΘU

for all s ∈ [d]. The proof of Proposition 21 gives F ∥Ds,q ∥L2 ≤ Cℓ A1/2 q ,

and also gives

F ∥Ds,q ∥L2 ≤ Cℓ A1/2 q ΘU ,

Using (48), we have

s ∈ [r], s ∈ [d] \ [r].

A1/2 ≤ ∥f ∗ ∥L2 (γd ) . q

Therefore, we have (Γq )s,t − (Ξq )s,t ≤ Cℓ ∥f ∗ ∥2L2 (γd ) ΘU 1{s≤r or t≤r} + Θ2U 1{s>r, t>r} . Combining this estimate with Proposition 21 gives e q )s,t ≤ Cℓ ∥f ∗ ∥2 (Γq )s,t − (G ΘU 1{s≤r or t≤r} + Θ2U 1{s>r, t>r} . L2 (γd )

If ΘU is sufficiently small, then Lemma 20 allows ∥f ∗ ∥2L2 (γd ) to be replaced by ∥f ∗ ∥2L2 (τd ) . The uniform version with arbitrary inactive unit vectors follows from the same unit-vector derivative estimate and the uniform inactive-direction estimate in Proposition 19. 53

I.10

Proof of Theorem 23

Since P≤L =

L X

Pq ,

q=0

the derivative (u[s] )⊤ ∇Pq is homogeneous multilinear of degree q − 1. Therefore the contributions from distinct degrees are orthogonal in L2 (τd ). Hence, for q ̸= q ′ , we have D E (u[s] )⊤ ∇Pq , (u[t] )⊤ ∇Pq′ = 0. L2 (τd )

It follows that (Γ≤L )s,t =

L X

(Γq )s,t .

(77)

e q )s,t . (G

(78)

q=1

On the Gaussian side, (50) gives e ≤L )s,t = (G

L X q=1

Subtracting (78) from (77), applying Proposition 22 for each q = 1, . . . , L, and absorbing the finite number of degrees into Cℓ proves the entrywise bounds. The same summation gives the stated uniform bounds involving arbitrary inactive unit vectors. It remains to pass from the entrywise comparison to the block-operator comparison. In the orthonormal basis {u[1] , . . . , u[d] }, e ≤L . Hence the desired operator the matrix M≤L is represented by Γ≤L , while U ⊤ ΣL U is represented by G norm equals e ≤L ∥op . ∥Γ≤L − G The active-active block has size r×r. Its entrywise bound therefore contributes at most Cℓ r ∥f ∗ ∥2L2 (τd ) ΘU to the operator norm. For the active-inactive block, the uniform inactive-direction estimate implies that, for every unit vector v ∈ U ⊥ and every s ∈ [r], the corresponding mixed bilinear form is bounded by Cℓ ∥f ∗ ∥2L2 (τd ) ΘU . Therefore the active-inactive block has operator norm at most √ Cℓ r ∥f ∗ ∥2L2 (τd ) ΘU . This is absorbed by Cℓ r ∥f ∗ ∥2L2 (τd ) ΘU . For the inactive-inactive block, the uniform inactive-inactive estimate gives operator norm at most Cℓ ∥f ∗ ∥2L2 (τd ) Θ2U . Since ΘU ≤ 1 in the regime of interest and r ≥ 1, this is also absorbed by Cℓ ∥f ∗ ∥2L2 (τd ) rΘU . Combining the three block estimates gives M≤L − U ⊤ ΣL U op ≤ Cℓ ∥f ∗ ∥2L2 (τd ) rΘU . Since ΘU = rΘU , this proves (51).

Proof of Theorem 25

Extend u[1] , . . . , u[r] to an orthonormal basis u[1] , . . . , u[d] of Rd , with u[r+1] . . . , u[d] ⊆ U ⊥ . For any u[s] , u[t] ∈ {u[i] }di=1 , we focus on the estimate n d h i 1 X X [s] ˆ (k) [t] (u[s] )⊤ En ∇fˆ(∇fˆ)⊤ u[t] = ui ∂i f (x )∂j fˆ(x(k) )uj nd2 i,j=1 k=1

1 ⊤ −1 y Kλ Diag(Xu[s] )(K ′ )2 Diag(Xu[t] )Kλ−1 y. = nd2 By the assumption H(f ∗ ) ∈ L2 (τd ), we can expand H(f ∗ ) into Fourier-Walsh basis: X H(fU∗ )(x(i) ) = ϕ(x(i) )⊤ c = cS ϕS (x(i) ),

(79)

(80)

S⊆[d],|S|≤ℓ

for some coefficients cS := Ex∼τd H(fU∗ )xS ,

(81)

S ⊆ [d],

which satisfies ∥c∥2 < ∞. Since E[yi2 ] = E|H(fU∗ )|2 + σε2 we deduce the estimate 2

E[∥y∥22 ] = n(E|H(fU∗ )|2 + σε2 ) = n(∥fU∗ ∥L2 (τd ) + σε2 ). Markov’s inequality subsequently shows 2

∥y∥ = Od,P (n log(n) · (∥fU∗ ∥L2 (τd ) + σε2 )). We summarize this observation in the following proposition. 2

Proposition 26. With probability at least 1 − 1/ log(n) we have ∥y∥ ≤ n log(n)(∥fU∗ ∥L2 (τd ) + σε2 ). We omit the subscript U in fU∗ for simplicity. Applying Lemma 12 with the two kernels K and K ′ shows that for any ϵ > 0 there exist matrices ∆, ∆1 satisfying ∥∆∥op , ∥∆1 ∥op = Od,P (d(δ−1)/2+ϵ ) and (82)

K = Φ≤p DΦ⊤ ≤p + ρIn + ∆1 , K

′

(83)

′

′ = Φ≤p D Φ⊤ ≤p + ρ In + ∆,

where we define ρ := g(1) − gp (1) and ρ′ := g ′ (1) − gp+1 (1), and the diagonal matrices D and D′ satisfy ∥DSk − g (k) (0)d−k I|Sk | ∥op = Od (d−k−1 ) and ∥DS′ k − g (k+1) (0)d−k I|Sk | ∥op = Od (d−k−1 ) for k = 0, . . . , p. With the expression (83) in place of K ′ , equation (79) becomes 2 1 ⊤ −1 ′ y Kλ Diag(Xu[s] ) Φ≤p D′ Φ⊤ Diag(Xu[t] )Kλ−1 y. ≤p + ρ In + ∆ 2 nd Letting β = Kλ−1 y and expanding the square gives h i (u[s] )⊤ En ∇fˆ(∇fˆ)⊤ u[t] =

2 1 ⊤ 1 β Diag(Xu[s] ) Φ≤p D′ Φ⊤ Diag(Xu[t] )β + 2 β ⊤ Diag(Xu[s] )(ρ′ In + ∆)2 Diag(Xu[t] )β ≤p 2 {z } nd | {z } nd | T2 (s,t)

T1 (s,t)

1 + 2 β ⊤ Diag(Xu[s] ) nd |

⊤

Φ≤p D′ Φ≤p (ρ′ In + ∆) + (ρ′ In + ∆) Φ≤p D′ Φ⊤ ≤p {z T3 (s,t)

Diag(Xu[t] )β }

The following claim is the centralhingredient iof the proof. Once T1 (s, t), T2 (s, t) and T3 (s, t) are estimated, the desired bound for (u[s] )⊤ En ∇fˆ(∇fˆ)⊤ u[t] follows immediately. Accordingly, the remainder of this section is devoted to the proof of this claim. In particular, the bound for T3 (s, t) follows immediately from the bound for T1 (s, t) and T2 (s, t) using Cauchy-Schwarz. Recall that (Γq )s,t := E (u[s] )⊤ ∇H(f ∗ )q (∇H(f ∗ )q )⊤ u[t] , s, t ∈ [d]. (84) To simplify the notation, for 0 ≤ q ≤ p + 1, define q q Aq (s, t) := (Γq )s,s + (Γq )t,t , and define A≤p (s, t) :=

Aq (s, t),

B≤p (s, t) :=

0≤q≤p

(Γq )s,s + (Γq )t,t .

0≤q≤p

Claim 2. The following estimates hold for any ϵ > 0 uniformly over all s, t ∈ [d]: 2 (a) nd1 2 |T2 (s, t)| = Od,P d−2+ϵ · ∥f ∗ ∥L2 + σε2 . (b)

1 nd2 T1 (s, t) −

− δ2 +ϵ

= Od,P d

0≤q≤p "

(Γq )s,t + C02 d2δ−2 (Γp+1 )s,t #

B≤p (s, t) + d2δ−2 ((Γp+1 )s,s + (Γp+1 )t,t )

1/2 1+δ ϵ 2 + Od,P d− 2 + 2 A≤p (s, t) ∥f ∗ ∥L2 + σε2 −4+3δ ϵ 1/2 2 + Od,P d 2 + 2 Ap+1 (s, t) ∥f ∗ ∥L2 + σε2 2 + Od,P d−1−δ+ϵ + d−2+δ+ϵ ∥f ∗ ∥L2 + σε2 , " # (c) nd1 2 |T3 (s, t)| = Od,P d−δ+ϵ B≤p (s, t) + d2δ−2 ((Γp+1 )s,s + (Γp+1 )t,t ) 1/2 1+3δ ϵ 2 + Od,P d− 2 + 2 A≤p (s, t) ∥f ∗ ∥L2 + σε2 −4+δ ϵ 1/2 2 + Od,P d 2 + 2 Ap+1 (s, t) ∥f ∗ ∥L2 + σε2 2 + Od,P d−1−2δ+ϵ + d−2+δ+ϵ ∥f ∗ ∥L2 + σε2 ,

J.1

Proof of Item (a) of Claim 2

Sub-multiplicativity of the operator norm directly implies 2

T2 (s, t) ≤ ∥ρ′ I + ∆∥op max |(u[s] )⊤ x(i) | max |(u[t] )⊤ x(i) | Kλ−1 y 2 . i

Note that the expression (82) directly implies Kλ−1 op ≤ 1/(ρ + λ − ∥∆1 ∥op ) = Od,P (1).

(85)

where we use the fact that D is positive semi-deifnite since D is diagonally dominant. Therefore, taking into account Proposition 26, we deduce ! r 2 −1 −1 ∗ 2 n log(n) ∥f ∥L2 + σε (86) Kλ y 2 ≤ Kλ op ∥y∥2 ≤ Od,P 56

By Hoeffding’s inequality we have for any t > 0, with probability at least 1 − 2e−t /2 over the random data x(i) , the following holds for any unit vector u: |u⊤ x(i) | ≤ t.

(87)

Letting t = log(d) and taking the union bound over {x(i) }ni=1 and {u[j] }dj=1 , we have max i∈[n],j∈[d] 2

(u[j] )⊤ x(i) = Od,P (log(d)).

With ∥ρ′ I + ∆∥op ≤ 2|ρ′ |2 + 2 ∥∆∥op = Od,P (1), consequently, we conclude 2 log (d) ∗ 2 1 2 2 −2+ϵ ∗ 2 2 |T (s, t)| = O · ∥f ∥ + σ = O d · ∥f ∥ + σ 2 d,P d,P ε ε L2 L2 nd2 d2

(88)

for any ϵ > 0 as desired.

J.2

Proof of Item (b) of Claim 2

Expanding the square and invoking Theorem 8, the term nd1 2 T1 (s, t) can be reduced to 1 1 [t] e ′ Φ⊤ T1 (s, t) = 2 β ⊤ Diag(Xu[s] )Φ≤p D′ (I + ∆)D ≤p Diag(Xu )β, nd2 d e where ∆

(89)

= Od,P (d−δ/2+ϵ ). To simplify the notation, letting [j] h[j] := d1 D′ Φ⊤ ≤p Diag(Xu )β

we write 1 e [t] . T1 (s, t) = (h[s] )⊤ h[t] + (h[s] )⊤ ∆h nd2

(90)

P d Observe that u X Φ≤p lies in the subspace span{Φ1 , Φ2 , . . . , Φp+1 }. Collecting coefficients with i i i=1 respect to this basis, we obtain ! d 1 ⊤ X β ui Xi Φ≤p D′ = β ⊤ Φ≤p+1 E, (91) d i=1 where the sparse matrix E ∈ RL(p+1)×L(p) corresponds to coefficients. A direct computation yields ( 1 ′ if S = J ⊔ {i} or J \ {i}; D · ui , ES,J = d J,J 0, otherwise. Here EJ\{i},J is induced by coordinate collapse. We observe that each column E:,J is a sparse vector with at most d non-zero entries. e [t] in equation (91) to bound T1 . Combining these two Next, we analyze (h[s] )⊤ h[t] and (h[s] )⊤ ∆h bounds completes the estimate for T1 . Specifically, we will prove the following holds with constant C0 = (ρ + λ)−1 g (p+1) (0):   X (h[s] )⊤ h[t] −  (Γq )s,t + C02 d2δ−2 (Γp+1 )s,t  0≤q≤p

= Od,P d

−δ+ 2ϵ

B≤p (s, t) −3+4δ 1/2 1+δ ϵ 2 + Od,P d− 2 + 2 A≤p (s, t) + d 2 Ap+1 (s, t) ∥f ∗ ∥L2 + σε2 2 + Od,P d−1−δ+ϵ + d−2+δ+ϵ ∥f ∗ ∥L2 + σε2

(92)

and e [t] = Od,P d− δ2 +ϵ B≤p (s, t) (h[s] )⊤ ∆h 3δ + Od,P d 2 −2+ϵ ((Γp+1 )s,s + (Γp+1 )t,t ) 3δ δ 2 + Od,P d−1− 2 +ϵ + d−2+ 2 +ϵ ∥f ∗ ∥L2 + σε2 .

(93)

Term (h[s] )⊤ h[t] . We decompose indices of vector (91) into j = 0, . . . , p − 1 and j = p, which correspondingly decomposes (h[s] )⊤ h[t] into two terms: D E D E [s] [t] [s] [t] (h[s] )⊤ h[t] = h≤p−1 , h≤p−1 + hSp , hSp . (94) Following the same convention, we let c≤q denote the vector [cS ]∪S:|S|≤q . To simplify the notation, corresponding to each direction j ∈ [d], we denote −1 r[j] := (E [j] )⊤ D≤p+1 c≤p+1 .

(95)

Next, we present two key lemmas that show h can be approximated by r. The proof appears in Appendix K.1 and K.2 respectively. Lemma 27. The following estimate holds for any ϵ > 0 and uniformly for any j ∈ [d] [j]

[j]

h≤p−1 − r≤p−1

[j]

= Od,P (d−2δ+ϵ ) r≤p−1

2 2

+ Od,P (d−1−δ+ϵ )(∥f ∗ ∥L2 + σε2 ).

(96)

Lemma 28. The following estimate holds for any ϵ > 0 and uniformly for any j ∈ [d] with constant C0 = (ρ + λ)−1 g (p+1) (0): [j]

[j]

= Od,P (d−2+ϵ ) rSp + Od,P (d−2+δ+ϵ )(∥f ∗ ∥L2 + σε2 ). 2 2 D E [s] [t] The following proposition connects rSq , rSq to Γs,t . The proof appears in Appendix K.3. hSp − C0 · dδ−1 rSp

(97)

Proposition 29. The following estimate holds for any q = 0, 1, . . . , p D E [s] [t] 2 rSq , rSq − (Γq+1 )s,t = Od (d−2 )Aq+1 (s, t) ∥f ∗ ∥L2 + Od (d−4 ) ∥f ∗ ∥L2 ,

(98)

Then, we have the following estimate. The proof appears in Appendix K.4. Proposition 30. for any ϵ > 0, the following estimates hold: D E X [t] [s] h≤p−1 , h≤p−1 − (Γq )s,t 0≤q≤p −δ+ 2ϵ

= Od,P d

B≤p (s, t) 1+δ ϵ 1/2 2 + Od,P d− 2 + 2 A≤p (s, t) ∥f ∗ ∥L2 + σε2 2 + Od,P d−1−δ+ϵ ∥f ∗ ∥L2 + σε2 ,

(99)

D E [s] [t] hSp , hSp − C02 d2δ−2 (Γp+1 )s,t −4+3δ ϵ 1/2 2 = Od,P d 2 + 2 Ap+1 (s, t) ∥f ∗ ∥L2 + σε2 2 + Od,P d−2+δ+ϵ ∥f ∗ ∥L2 + σε2 .

(100)

and

Plugging the above two estimates into equation (94) finishes the bound for (92). 58

e [t] . Term (h[s] )⊤ ∆h

First, applying (94) with s = t, together with Proposition 30, gives 2

1 − B≤p (s, s) − C02 d2δ−2 (Γp+1 )s,s 2 −δ+ 2ϵ ≤ Od,P d B≤p (s, s) 1/2 1+δ ϵ 2 + Od,P d− 2 + 2 A≤p (s, s) ∥f ∗ ∥L2 + σε2 −4+3δ ϵ 1/2 2 + Od,P d 2 + 2 Ap+1 (s, s) ∥f ∗ ∥L2 + σε2 2 + Od,P d−1−δ+ϵ + d−2+δ+ϵ ∥f ∗ ∥L2 + σε2 . h[s]

(101)

P Here the factor 1/2 appears because B≤p (s, s) = 2 0≤q≤p (Γq )s,s . Next, since 2  X q (Γq )s,s  ≤ 2(p + 1)B≤p (s, s), A≤p (s, s)2 = 4  0≤q≤p

the weighted AM-GM inequality implies, after applying the estimate with a smaller ϵ and then renaming it as ϵ, that 1+δ ϵ 1/2 2 Od,P d− 2 + 2 A≤p (s, s) ∥f ∗ ∥L2 + σε2 2 = od,P (1)B≤p (s, s) + Od,P d−1−δ+ϵ ∥f ∗ ∥L2 + σε2 . (102) Similarly, since Ap+1 (s, s) = 2

p (Γp+1 )s,s , another application of weighted AM-GM gives

−4+3δ ϵ 1/2 2 Od,P d 2 + 2 Ap+1 (s, s) ∥f ∗ ∥L2 + σε2 2 = od,P (1)d2δ−2 (Γp+1 )s,s + Od,P d−2+δ+ϵ ∥f ∗ ∥L2 + σε2 .

(103)

Therefore, substituting (102) and (103) into (101), we obtain 2 1 + od,P (1) B≤p (s, s) + C02 + od,P (1) d2δ−2 (Γp+1 )s,s h[s] = 2 2 2 + Od,P d−1−δ+ϵ + d−2+δ+ϵ ∥f ∗ ∥L2 + σε2 .

(104)

Likewise, replacing s by t, we have 2 1 h[t] = + od,P (1) B≤p (t, t) + C02 + od,P (1) d2δ−2 (Γp+1 )t,t 2 2 2 + Od,P d−1−δ+ϵ + d−2+δ+ϵ ∥f ∗ ∥L2 + σε2 .

(105)

Consequently, adding (104) and (105), and using 12 B≤p (s, s) + 12 B≤p (t, t) = B≤p (s, t), we get h[s]

+ h[t]

2 2

= (1 + od,P (1)) B≤p (s, t) + C02 + od,P (1) d2δ−2 ((Γp+1 )s,s + (Γp+1 )t,t ) 2 + Od,P d−1−δ+ϵ + d−2+δ+ϵ ∥f ∗ ∥L2 + σε2 .

(106)

Finally, by Cauchy’s inequality and AM-GM, the following holds [s] ⊤ e

(h ) ∆h

[t]

e ≤ ∆

[s]

[t]

h op

h 2

1 e ≤ ∆ 2 op 2

[s]

[t]

+ h 2

e Thus, using ∆

δ = Od,P d− 2 +ϵ , and substituting (106), we obtain e [t] = Od,P d− δ2 +ϵ B≤p (s, t) (h[s] )⊤ ∆h 3δ + Od,P d 2 −2+ϵ ((Γp+1 )s,s + (Γp+1 )t,t ) 3δ δ 2 + Od,P d−1− 2 +ϵ + d−2+ 2 +ϵ ∥f ∗ ∥L2 + σε2 .

(107)

which gives (93).

Proof of lemmas and propositions for the main results

K.1

Proof of Lemma 27

[t] For notational simplicity, we suppress Sp the superscript [t] in u . A direct coordinate calculation gives the following identity: for every J ∈ j=0 Sj ,

 ui , if S = J ⊔ {i},    1 |S| −1 −1 |S|+2 (D≤p+1 E:,J )S = g (0) g (0) ui , if S = J \ {i}, 2  d   0, otherwise.

(108)

Consequently, this identity implies the uniform bound −1 D≤p+1 E:,J

We next restrict to J ∈ of h satisfies

(109)

= Od (1).

Sp−1

j=0 Sj . In this case, ESp+1 ,J = 0. Therefore, the corresponding coordinate

hJ = β ⊤ Φ≤p+1 E:,J = β ⊤ Φ≤p E:L(p),J + ΦSp+1 ESp+1 ,J

= β ⊤ Φ≤p E:L(p),J . Since D≤p+1 is diagonal, we define −1 −1 êJ := (D≤p+1 E:,J ):L(p) = D≤p E:L(p),J ∈ RL(p) .

(110)

Thus, when the context is clear, we write D in place of D≤p . By Lemma 13, applied with m = 2p + 1, the regularized kernel admits the refined decomposition (111)

Kλ = Φ≤p DΦ⊤ ≤p + (ρp,2p+1 + λ)H, where H is given by  H = In + (ρp,2p+1 + λ)−1 

2p+1 X

 θk off ΦSk ΦSk + ∆ .

(112)

(113)

⊤

k=p+1

Moreover, the perturbation ∆ satisfies ∥∆∥op = Od,P

log n d

2p+2−p−δ −ϵ 2

= Od,P (n−1/2 ).

Here the coefficients satisfy θk = g (k) (0)d−k + Od (d−k−1 ) for k = p + 1, . . . , 2p + 1, and ρp,2p+1 = g(1) − gp (1) + Od (d−1 ). For notational simplicity, write ρ := ρp,2p+1 . Then H can be written as (114)

H = In + ∆1 , where the perturbation is  ∆1 := (ρ + λ)−1 

2p+1 X

 (115)

⊤

θk off ΦSk ΦSk + ∆ .

k=p+1

By Lemma 12, this perturbation satisfies ∥∆1 ∥op = Od,P d(δ−1)/2+ϵ .

(116)

Λ := (nD)−1 .

(117)

Set

Since E:L(p),J = DêJ , the Sherman–Morrison–Woodbury formula gives hJ = β ⊤ Φ≤p DêJ = (f ∗ (X) + ε)⊤ Kλ−1 Φ≤p DêJ  =

−1 −1 Φ⊤ Φ≤p  ≤p H

1 ∗  (f (X) + ε)⊤ H −1 Φ≤p (ρ + λ)Λ + n | {z

=:G

 }

êJ

1 1 = ε⊤ H −1 Φ≤p G−1 êJ + f ∗ (X)⊤ H −1 Φ≤p G−1 êJ . n n

(118)

We will use the following norm bounds. Their proof is given in Appendix K.5. Claim 3. The following bounds hold: H −1 op = Od,P (1),

(119)

G−1 op = Od,P (1).

We first bound the noise term in (118). Define (3)

hJ :=

1 ⊤ −1 ε H Φ≤p G−1 êJ . n

(120)

Conditioning on X and using the independence of ε, we obtain 2 σ 2 2 (3) Eε hJ = ε2 H −1 Φ≤p G−1 êJ 2 n σε2 2 2 2 2 ≤ 2 H −1 op ∥Φ≤p ∥op G−1 op ∥êJ ∥2 n = Od,P (n−1 ) σε2 . In the last line, we used (119), Lemma 7, and (109). More explicitly, these bounds give H −1 op = Od,P (1),

G−1 op = Od,P (1), 61

∥Φ≤p ∥op = Od,P (n),

∥êJ ∥2 = Od (1).

After summing over all J ∈ Sj , j = 0, . . . , p − 1, and applying Markov’s inequality, we get, for every ϵ > 0, p−1 X 2 X (3) hJ = Od,P d−1−δ+ϵ σε2 .

(121)

j=0 J∈Sj

We next decompose the signal part according to degree. Namely, write (122)

f ∗ (X) = Φ≤p c≤p + Φ>p c>p . Substituting this decomposition into the second term in (118) gives 1 ∗ 1 ⊤ −1 f (X)⊤ H −1 Φ≤p G−1 êJ = c⊤ Φ≤p G−1 êJ ≤p Φ≤p H n n | {z } (1)

=:hJ

1 ⊤ ⊤ −1 c>p Φ>p H Φ≤p G−1 êJ . n {z } | (2)

=:hJ

The remaining estimates are supplied by the following two propositions. Their proofs are given in Appendices K.6 and K.7, respectively. Proposition 31. For every ϵ > 0, the following estimate holds uniformly over {u[j] }dj=1 : X

(1)

hJ − c⊤ ≤p êJ

2 2 = Od,P d−2δ+ϵ ∥r≤p−1 ∥2 + Od,P d−p−δ+ϵ ∥f ∗ ∥L2 .

(123)

J:|J|≤p−1

Proposition 32. For every ϵ > 0, the following estimate holds uniformly over {u[j] }dj=1 :

(2)

2 = Od,P d−1−δ+ϵ ∥f ∗ ∥L2 .

(124)

J:|J|≤p−1 (1)

(2)

(3)

We now combine the three pieces. Since hJ = hJ + hJ + hJ for |J| ≤ p − 1, and since rJ = c⊤ ≤p êJ , the elementary inequality (a + b + c)2 ≤ 3a2 + 3b2 + 3c2 gives 2 ∥h≤p−1 − r≤p−1 ∥2 =

p−1 X X

hJ − c⊤ ≤p êJ

j=0 J∈Sj

≤3

(1)

hJ − c⊤ ≤p êJ

J:|J|≤p−1

(2)

J:|J|≤p−1

(3)

J:|J|≤p−1

Using (123), (124), and (121), we therefore obtain 2 2 2 ∥h≤p−1 − r≤p−1 ∥2 = Od,P d−2δ+ϵ ∥r≤p−1 ∥2 + Od,P d−p−δ+ϵ ∥f ∗ ∥L2 2 + Od,P d−1−δ+ϵ ∥f ∗ ∥L2 + σε2 . 2

Since p ≥ 1, the term d−p−δ+ϵ ∥f ∗ ∥L2 is absorbed by the d−1−δ+ϵ ∥f ∗ ∥L2 term. Hence 2 2 2 ∥h≤p−1 − r≤p−1 ∥2 = Od,P d−2δ+ϵ ∥r≤p−1 ∥2 + Od,P d−1−δ+ϵ ∥f ∗ ∥L2 + σε2 . This proves the desired estimate. 62

K.2

Proof of Lemma 28

Fix J ∈ Sp . We begin by substituting the decomposition of y into the definition of hJ . Namely, we use y = ΦSp+1 cSp+1 + ΦQ cQ + ε. Since

hJ = β ⊤ ΦSp+1 ESp+1 ,J ,

the substitution gives −1 ⊤ hJ = c⊤ Sp+1 ΦSp+1 Kλ ΦSp+1 ESp+1 ,J ⊤ −1 ⊤ −1 + c⊤ Q ΦQ Kλ ΦSp+1 ESp+1 ,J + ε Kλ ΦSp+1 ESp+1 ,J .

(125)

To simplify notation, set eJ := ESp+1 ,J . We first estimate the noise contribution. Taking expectation over ε, conditionally on the design, gives h 2 i 2 Eε ε⊤ Kλ−1 ΦSp+1 eJ = σε2 Kλ−1 ΦSp+1 eJ 2 . Using the estimates ∥Kλ−1 ∥op = Od,P (1) from (85) and ∥ΦSp+1 eJ ∥22 = Od,P (dp+δ+ϵ )∥eJ ∥22 from (26), we obtain h 2 i Eε ε⊤ Kλ−1 ΦSp+1 eJ = Od,P d−p−2+δ+ϵ σε2 . After summing over J ∈ Sp and applying Markov’s inequality, the preceding bound yields X 2 ε⊤ Kλ−1 ΦSp+1 eJ = Od,P d−2+δ+ϵ σε2 .

(126)

J∈Sp

We next estimate the first deterministic term in (125). The proof of the following estimate is given in Appendix K.8. Proposition 33. For every ϵ > 0, the following estimate holds: 2 X 1 1 ⊤ −1 ⊤ c⊤ c Φ e K Φ e − J Sp+1 J n Sp+1 Sp+1 λ ρ + λ Sp+1 J∈Sp 2 X ⊤ 2 = Od,P d−2δ+ϵ cSp+1 eJ + Od d−2p−2−δ+ϵ cSp+1 2 .

(127)

J∈Sp

We then estimate the second deterministic term in (125). The proof of the next estimate is given in Appendix K.9. Proposition 34. For every ϵ > 0, the following estimate holds: 2 X 1 2 ⊤ −1 c⊤ Φ K Φ e = Od d−2p−2−δ+ϵ ∥cQ ∥2 . Sp+1 J Q Q λ n

(128)

J∈Sp

We now combine the three estimates. First, applying (a + b + c)2 ≤ 3(a2 + b2 + c2 ) to the decomposition (125) gives 2 X n ⊤ hJ − c eJ ρ + λ Sp+1 J∈Sp 2 X X 2 n ⊤ −1 ⊤ ⊤ −1 ⊤ ⊤ ≤3 cQ ΦQ Kλ ΦSp+1 eJ + 3 cSp+1 ΦSp+1 Kλ ΦSp+1 eJ − c eJ ρ + λ Sp+1 J∈Sp J∈Sp X 2 +3 ε⊤ Kλ−1 ΦSp+1 eJ . (129) J∈Sp

Multiplying the estimates in Propositions 33 and 34 by n2 , and then using (126), gives 2 X n ⊤ hJ − c eJ ρ + λ Sp+1 J∈Sp 2 X ⊤ = Od,P d−2δ+ϵ n cSp+1 eJ + Od d−2+δ+ϵ ∥f ∗ ∥2L2 + σε2 .

(130)

J∈Sp

It remains to rewrite the leading deterministic term in the desired normalization. For every J ∈ Sp , the diagonal expansion gives DJ,J = g (p+1) (0)d−p−1 + Od (d−p−2 ). Therefore, multiplying by n and using the scaling of n, we have nDJ,J = g (p+1) (0)d−1+δ + Od (d−2+δ ). Consequently, the following identity holds: −1 ⊤ n c⊤ Sp+1 eJ = nDJ,J DJ,J cSp+1 eJ . Substituting the preceding estimate for nDJ,J , we obtain −1 (p+1) n c⊤ (0)d−1+δ + Od (d−2+δ ) c⊤ Sp+1 eJ = g Sp+1 DSp+1 eJ . We now substitute this expression into (130). The result is 2 g (p+1) (0) −1+δ ⊤ −1 hJ − d cSp+1 DSp+1 eJ ρ+λ J∈Sp 2 X ⊤ = Od,P d−2+ϵ cSp+1 DS−1 e + Od d−2+δ+ϵ ∥f ∗ ∥2L2 + σε2 . J p+1 X

(131)

J∈Sp −1 Finally, define C0 := (ρ + λ)−1 g (p+1) (0). Also recall the notation rJ := c⊤ Sp+1 DSp+1 eJ . With these definitions, (131) becomes X 2 X 2 hJ − C0 d−1+δ rJ = Od,P d−2+ϵ rJ + Od d−2+δ+ϵ ∥f ∗ ∥2L2 + σε2 . J∈Sp

J∈Sp

This is the desired estimate.

K.3

Proof of Proposition 29

We set

P := H(fU∗ ),

Pq := [P ]deg q .

W also write Pq = H(fU∗ )q By definition, we have

ui cJ⊔{i} xJ = u⊤ ∇Pq+1 .

J:|J|=q i∈[d]\J

With this expression, we can deduce 

 h

E (u[s] )⊤ ∇Pq+1 (∇Pq+1 )⊤ u[t] =

X 

J:|J|=q

i∈[d]\J

[s]

ui cJ⊔{i}  

 X

i∈[d]\J

[t]

ui cJ⊔{i}  .

(132)

For any J ∈ Sq , we can write rJ as ⊤ rJ = E:,J D≤p+1 c≤p+1 X 1X ′ −1 DJ,J DJ\{i},J\{i} ui cJ\{i} , = ui cJ⊔{i} + d i∈J

i∈[d]\J

which implies rJ −

X 1X ′ −1 ui cJ\{i} |. DJ,J DJ\{i},J\{i} ui cJ\{i} = Od (d−2 )| d

ui cJ⊔{i} =

(133)

i∈J

i∈[d]\J

Using basic inequality for any a, b, c, d ∈ Rn (134)

| ⟨a, b⟩ − ⟨c, d⟩ | ≤ ∥d∥2 ∥a − c∥2 + ∥c∥2 ∥b − d∥2 + ∥a − c∥2 ∥b − d∥2 , Eq. (132) and (133) together gives D E [s] [t] rSq , rSq − (Γq+1 )s,t =

h i [s] [t] rJ rJ − E (u[s] )⊤ ∇Pq+1 (∇Pq+1 )⊤ u[t]

X J:|J|=q

v 2 v  u u u X u X X u u [s]  ≤ Od (d−2 )t ui cJ⊔{i}  t J:|J|=q

i∈[d]\J

!2 X

[s] ui cJ\{i}

i∈J

J:|J|=q

v 2 v  u u u X u X X u u [t]  ui cJ⊔{i}  t + Od (d−2 )t

i∈[d]\J

i∈J

J:|J|=q

v u u X −4 u + Od (d )t J:|J|=q

J:|J|=q

v !2 u u X X [s] u ui cJ\{i} t

i∈J

J:|J|=q

[t] ui cJ\{i}

!2 [t] ui cJ\{i}

By Cauchy-Schwarz inequality, we have v u !2 s u X X X X u 2 t cJ\{i} ≤ ∥f ∗ ∥L2 . ui cJ\{i} ≤ i∈J

J:|J|=q

(135)

(136)

J:|J|=q i∈J

Plugging this back into equation (135) yields D

E [s] [t] rSq , rSq − (Γq+1 )s,t =

h i [s] [t] rJ rJ − E (u[s] )⊤ ∇Pq+1 (∇Pq+1 )⊤ u[t]

J:|J|=q

≤ Od (d

−2

q q 2 ) (Γq+1 )s,s + (Γq+1 )t,t ∥f ∗ ∥L2 + Od (d−4 ) ∥f ∗ ∥L2 ,

which completes the proof.

K.4

Proof of Proposition 30

We first record two residual-size estimates which follow from the diagonal case of Proposition 29. Indeed, by applying Proposition 29 with s = t, and then summing over the orthogonal components, we obtain [s]

r≤p−1

2 2

[t]

+ r≤p−1

= B≤p (s, t) + Od (d

2 2 −2

)A≤p (s, t) ∥f ∗ ∥L2 + Od (d−4 ) ∥f ∗ ∥L2 . 65

(137)

Moreover, taking square roots componentwise in the same diagonal estimate gives [t] [s] r≤p−1 + r≤p−1 = Od A≤p (s, t) + d−2 ∥f ∗ ∥L2 .

(138)

Similarly, for the top component, the same argument gives [t]

[s]

rSp

+ rSp

(139)

= Od Ap+1 (s, t) + d−2 ∥f ∗ ∥L2 .

We now prove (99). First, for u ∈ {s, t}, the approximation estimate (96) gives 1/2 1+δ ϵ ϵ [u] [u] [u] 2 h≤p−1 − r≤p−1 = Od,P d−δ+ 2 r≤p−1 + Od,P d− 2 + 2 ∥f ∗ ∥L2 + σε2 . 2

Therefore, applying (134) with [s]

[t]

a = h≤p−1 , we obtain

[s]

b = h≤p−1 ,

[t]

c = r≤p−1 ,

d = r≤p−1 ,

D E D E [s] [t] [s] [t] h≤p−1 , h≤p−1 − r≤p−1 , r≤p−1 [t]

[s] [s] [s] [t] [t] h≤p−1 − r≤p−1 + r≤p−1 h≤p−1 − r≤p−1 2 2 2 2 [s] [s] [t] [t] h≤p−1 − r≤p−1 h≤p−1 − r≤p−1 . 2 2

≤ r≤p−1 +

Consequently, substituting the preceding approximation estimate and using AM-GM yields D E D E [s] [t] [s] [t] h≤p−1 , h≤p−1 − r≤p−1 , r≤p−1 2 2 [s] [t] −δ+ 2ϵ r≤p−1 + r≤p−1 ≤ Od,P d 2

1+δ

[s]

[t]

+ Od,P d− 2 + 2 r≤p−1 + r≤p−1 2 2 −1−δ+ϵ 2 ∗ 2 + Od,P d ∥f ∥L2 + σε .

1/2 2 ∥f ∗ ∥L2 + σε2

Next, using (137) and (138), and using 1/2 2 ∥f ∗ ∥L2 ≤ ∥f ∗ ∥L2 + σε2 , we get

D E D E [s] [t] [s] [t] h≤p−1 , h≤p−1 − r≤p−1 , r≤p−1 ϵ = Od,P d−δ+ 2 B≤p (s, t) 1+δ ϵ 1/2 2 + Od,P d− 2 + 2 A≤p (s, t) ∥f ∗ ∥L2 + σε2 2 + Od,P d−1−δ+ϵ ∥f ∗ ∥L2 + σε2 . 2

Here the first term is kept as B≤p (s, t), rather than being replaced by A≤p (s, t)(∥f ∗ ∥L2 + σε2 ), because A≤p (s, t) already has the scale of ∥f ∗ ∥L2 . On the other hand, by the orthogonality of the components and Proposition 29, we also have D

E X [s] [t] r≤p−1 , r≤p−1 − (Γq )s,t 0≤q≤p

≤

E D [t] [s] rSq−1 , rSq−1 − (Γq )s,t

0≤q≤p 2

= Od (d−2 )A≤p (s, t) ∥f ∗ ∥L2 + Od (d−4 ) ∥f ∗ ∥L2 . 66

Since 0 < δ < 1, the last display may be weakened as D

E X [s] [t] r≤p−1 , r≤p−1 − (Γq )s,t 0≤q≤p

1/2 1+δ 2 = Od d− 2 + 2 A≤p (s, t) ∥f ∗ ∥L2 + σε2 2 + Od d−1−δ+ϵ ∥f ∗ ∥L2 + σε2 .

Finally, combining the previous two estimates by the triangle inequality gives (99). We next prove (100). First, for u ∈ {s, t}, the Sp -component approximation gives δ−2 ϵ 1/2 [u] [u] 2 hSp − C0 dδ−1 rSp = Od,P d 2 + 2 ∥f ∗ ∥L2 + σε2 . 2

Therefore, applying (134) with [s]

[t]

a = hSp , we obtain

[s]

c = C0 dδ−1 rSp ,

b = hSp ,

[t]

d = C0 dδ−1 rSp ,

E E D [s] [t] [s] [t] hSp , hSp − C02 d2δ−2 rSp , rSp [t]

≤ C0 dδ−1 rSp δ−1

+ C0 d +

[s]

hSp − C0 dδ−1 rSp

[s] rSp

2 [t] δ−1 [t] hSp − C0 d rSp

2 [s] δ−1 [s] hSp − C0 d rSp

2 [t] δ−1 [t] hSp − C0 d rSp

. 2

Consequently, substituting the preceding Sp -component approximation gives D E D E [s] [t] [s] [t] hSp , hSp − C02 d2δ−2 rSp , rSp 1/2 δ−2 ϵ [s] [t] 2 ≤ Od,P dδ−1 d 2 + 2 rSp + rSp ∥f ∗ ∥L2 + σε2 2 2 ∗ 2 −2+δ+ϵ 2 + Od,P d ∥f ∥L2 + σε . Next, using (139), we get E E D D [s] [t] [s] [t] hSp , hSp − C02 d2δ−2 rSp , rSp 1/2 −4+3δ ϵ 2 = Od,P d 2 + 2 Ap+1 (s, t) ∥f ∗ ∥L2 + σε2 2 + Od,P d−2+δ+ϵ ∥f ∗ ∥L2 + σε2 . D E [s] [t] It remains to replace rSp , rSp by (Γp+1 )s,t . By Proposition 29 with q = p, we have D E [s] [t] C02 d2δ−2 rSp , rSp − (Γp+1 )s,t 2 ≤ C02 d2δ−2 Od (d−2 )Ap+1 (s, t) ∥f ∗ ∥L2 + Od (d−4 ) ∥f ∗ ∥L2 . Since 0 < δ < 1, the preceding display implies D E [s] [t] C02 d2δ−2 rSp , rSp − (Γp+1 )s,t 1/2 −4+3δ 2 Ap+1 (s, t) ∥f ∗ ∥L2 + σε2 = Od d 2 2 + Od d−2+δ ∥f ∗ ∥L2 + σε2 . Finally, combining the last two estimates by the triangle inequality proves (100), and the proof is complete. 67

K.5

Proof of Claim 3

Proof. The following estimate for any ∆ with ∥∆∥op = od,P (1): (140)

(I + ∆)−1 − I op = Od,P (∥∆∥op ), which can be deduced using Neumann series. Based on the definition of D, we immediately have (nD)−1 op = Od (d−δ ). Notice that

Φ⊤ ≤p Φ≤p − In n

= op

Od,P (d−δ/2+ϵ ) (Theorem 8) and (141)

H −1 − In op = Od (∥∆1 ∥op ) = Od,P (d(δ−1)/2+ϵ ), which follows from Eq. (140). Thus we deduce −1 Φ⊤ Φ≤p ≤p H −I n

= op

≤

−1 Φ⊤ Φ⊤ − I)Φ≤p ≤p Φ≤p ≤p (H + −I n n

Φ⊤ ≤p Φ≤p

−I

= Od,P (d

Φ≤p √ n

+ op

−δ/2+ϵ

(δ−1)/2+ϵ

2 op

H −1 − I op (142)

With this bound, we can deduce ∥G − I∥op = Od (d−δ ) + Od,P (d−δ/2+ϵ + d(δ−1)/2+ϵ ) = Od,P (d−δ/2+ϵ + d(δ−1)/2+ϵ ). Invoking Eq. (140) again we obtain (143)

G−1 − I op = Od,P (d−δ/2+ϵ + d(δ−1)/2+ϵ ).

K.6

Proof of Proposition 31

First, we express H −1 and G−1 in (118) using feature matrices. The proof appears in Appendix K.10. Proposition 35. The following equalities hold:  q 2p+1 cH X X c j   +∆1 , H −1 = In + − j off ΦSj Φ⊤ Sj d q=1 j=p+1 | {z }

(144)

:=∆H

and −1

= I|L(p)| +

cG X

" −(ρ + λ)Λ + n1 off

Φ⊤ ≤p Φ≤p

r=1

Φ⊤ Φ≤p ≤p + √ ∆H √ n n

+∆2 , }

:=∆G

where cj , cG , cH > 0 are constants and matrices ∆1 , ∆2 satisfy 1

∥∆1 ∥op , ∥∆2 ∥op = Od,P (n− 2 ).

#r (145)

To simplify the notation, we write AH :=

1 ⊤ Φ ∆H Φ≤p , n ≤p

AG :=

1 ⊤ Φ Φ≤p ∆G , n ≤p

AHG :=

1 ⊤ Φ ∆H Φ≤p ∆G n ≤p

(146)

−1 Plugging the expressions of H −1 = I + ∆H + ∆1 and G−1 = I + ∆G + ∆2 into Φ⊤ Φ≤p G−1 and ≤p H expanding the powers, we see that it can be written as

1 −1 1 ⊤ Φ≤p G−1 = Φ⊤ Φ≤p + AH + AG + AHG n Φ≤p H n ≤p (147)

+ n1 (Φ≤p )⊤ H −1 (Φ≤p )∆2 + n1 (Φ≤p )⊤ ∆1 (Φ≤p )G−1 . | {z } :=∆3

We can immediately deduce that ∆3 satisfies 1

∥∆3 ∥op = Od,P (1)(∥∆2 ∥op + ∥∆4 ∥op ) = Od,P (n− 2 ),

(148)

√ where we use the estimate of H −1 op = Od,P (1) (141), G−1 op = Od,P (1) (143) and ∥Φ≤p / n∥op = Od,P (1) (Lemma 7) The following proposition estimates AH , AG and AHG . The proof appears in Appendix K.11. Proposition 36. Let PG (Λ) :=

cG X

r −(ρ + λ)Λ .

r=1

For all N ∈ {AH , AG − PG (Λ), AHG }, the following equalities hold: ! 2 2 −4−2δ E (c⊤ = Od (d−2δ )(c⊤ )· ≤p N êJ ) ≤p êJ ) + Od (d

+ Od,P (n−1 ) ∥c≤p ∥2 .

c2J\{i}

(149)

i∈J (1)

Then, we plug the decomposition (147) into the expression of hJ and using Cauchy-Schwartz inequality we obtain 2 (1) hJ − ⟨c≤p , (I + PG (Λ)) êJ ⟩ 2 2 1 ⊤ ≤ 5 c⊤ + 5 c⊤ ≤p n Φ≤p Φ≤p + PG (Λ) êJ − ⟨c≤p , (I + PG (Λ)) êJ ⟩ ≤p AH êJ 2 2 2 + 5 c⊤ + 5 c⊤ + 5 c⊤ . (150) ≤p AG êJ ≤p (AHG − PG (Λ)) (êJ ) ≤p ∆3 (êJ ) 2 Using equation (148), the last term c⊤ can be bounded by ≤p ∆3 (êJ ) 2

(151)

∥∆3 ∥op · ∥c≤p ∥2 · ∥êJ ∥2 = Od,P (n−1 ) ∥c≤p ∥2 . Taking the expectation of the first term and expanding the square gives h 2 i 1 ⊤ E c⊤ Φ Φ + P (Λ) ê − ⟨c , (I + P (Λ)) ê ⟩ G J ≤p G J ≤p n ≤p ≤p h i 2 2 1 ⊤ = E c⊤ − c⊤ ≤p (I + PG (Λ))êJ ≤p n Φ≤p Φ≤p + PG (Λ) êJ (a)

= c⊤ ≤p (I + PG (Λ))êJ

+ Od ( n1 ) ∥c≤p ∥2 ∥(I + PG (Λ))êJ ∥2 − c⊤ ≤p (I + PG (Λ))êJ

2 (152)

= Od ( n1 ) ∥c≤p ∥2 , 69

where (a) follows from Lemma 39. For the second, third and fourth terms in equation (150), we take the expectation and the bound (149) immediately applies. Applying Proposition 43 shows the bound (149) and (152) hold uniformly for all u[1] , . . . u[d] , with an extra multiplicative factor dϵ . Summing over J with |J| ≤ p − 1 and applying Markov’s inequality, we have with probability 1 − 1/ log(d), the following holds 2 X 1 ⊤ Φ Φ≤p + AH + AG + AHG − ⟨c≤p , (I + PG (Λ))êJ ⟩ n ≤p J:|J|≤p−1

= Od,P (d−1−δ+ϵ ) log(d) ∥c≤p ∥2 . Combining with equation (151), we conclude 2 X (1) 2 2 hJ − ⟨c≤p , (I + PG (Λ))êJ ⟩ = Od (d−1−δ+ϵ ) ∥c≤p ∥2 = Od,P (d−1−δ+ϵ ) ∥f ∗ ∥L2 , J:|J|≤p−1

where we absorb the logarithmic actor into dϵ . Lastly, using the Cauchy-Schwarz inequality, we obtain 2 X (1) X 2 2 hJ − ⟨c≤p , êJ ⟩ = 2 c⊤ + Od,P (d−1−δ+ϵ ) ∥f ∗ ∥L2 . ≤p PG (Λ)êJ J:|J|≤p−1

(153)

J:|J|≤p−1

2 The following claim shows that c⊤ is of a smaller order then combining it with the equa≤p PG (Λ)êJ tion (153) completes the proof. The proof of the claim appears in Appendix K.12. Claim 4. The following estimate holds: c⊤ ≤p PG (Λ)êJ

! 2

= Od (d−2δ ) (c≤p êJ ) + Od (d−8−2δ )

c2J\{i}

(154)

i∈J

K.7

Proof of Proposition 32

The proof follows analogously to the proof of Proposition 31. For (2)

hJ =

1 ⊤ ⊤ −1 c Φ H Φ≤p G−1 êJ , n >p >p

plugging the expression of H −1 (144) and G−1 (145) gives ηJ =

1 ⊤ ⊤ c Φ (I + ∆H + ∆1 )Φ≤p (I + ∆G + ∆2 )êJ n >p >p

As in the proof of Proposition 31 , the terms involving the remainder matrices ∆1 , ∆2 are easier and are absorbed into the same bound, so it suffices to treat the principal part XJ :=

1 ⊤ ⊤ c Φ (I + ∆H )Φ≤p (I + ∆G )êJ . n >p >p

We decompose both endpoint spaces by exact degree. On the low-degree side, we use the same notation as in the proof of Proposition 36: Wq := n−1/2 ΦSq ,

W = [W0 , . . . , Wp ],

Λ=

p X

β q Pq ,

βq :=

q=0

On the high-degree side, write Φ>p = [ΦSp+1 , . . . , ΦSℓ ],

Vj := |Sj |−1/2 ΦSj 70

(p + 1 ≤ j ≤ ℓ),

dq . n

and decompose c>p = (cSp+1 , . . . , cSℓ ). Exactly as in the proof of Proposition 36, every summand in ∆H is a finite product of matrices off(Vj Vj⊤ ) with deterministic coefficient Od (1), and every word in ∆G is obtained by expanding powers of B0 := −(ρ + λ)Λ,

B1 := off(W ⊤ W ),

B2 := W ⊤ ∆H W.

Hence, after expanding I + ∆H and I + ∆G , every summand contributing to XJ is of the form 1 ⊤ ⊤ c Φ H0 Φ≤p Ξ êJ , n Sj Sj where p + 1 ≤ j ≤ ℓ, H0 is either I or one of the words appearing in ∆H , and Ξ is either I or one of the words appearing in ∆G . We now expand the low-degree factor Φ≤p Ξ exactly as in the proof for the low-degree endpoint words, by inserting the exact-degree projectors Pq around every copy of Λ and using the decompositions of B1 and B2 . This yields a finite representation X XJ = αη (aSjη )⊤ Mη (Pqη êJ ) , (155) η∈Ω

where |Ω| = Op,cH ,cG (1), each Mη is a normalized chain whose families belong to S0 , . . . , Sℓ , and whose endpoint families are Sjη and Sqη . In particular, the endpoint families are always different, since jη > p whereas qη ≤ p. Moreover, each coefficient αη is deterministic and satisfies ! r |Sjη | . (156) |αη | = Od n Indeed, the outer factor n−1 Φ⊤ Sj ΦSq contributes 1 ⊤ Φ ΦS = n Sj q

|Sj | ⊤ V Wq , n j

while all remaining coefficients coming from ∆H are Od (1), and every additional factor coming from Λ is one of the βq ≤ 1. Fix one η ∈ Ω. Since the endpoint families of Mη are different, Lemma 39 gives h 2 i (1) (m+1) 2 −1 2 2 E c⊤ M (P ê ) ≤ O (nw w ) n ∥a ∥ ∥P ê ∥ (157) η qη J d Sjη 2 qη J 2 . Sjη Note that the weights satisfy w(1) = |Sjη |−1/2 , therefore we have

w(m+1) = n−1/2 ,

(nw(1) w(m+1) )2 =

n . |Sjη |

Substituting this into 157 yields h 2 i 1 2 2 E c⊤ M (P ∥c ê ) ≤ O ∥ ∥P ê ∥ η q J d S q J Sjη 2 . η jη 2 η |Sjη | Combining this with (156), we obtain h 2 i E αη (aSjη )⊤ Mη (Pqη êJ ) ≤ Od n−1 ∥aSjη ∥22 ∥Pqη êJ ∥22 . 71

(158)

Since |Ω| = Op,cH ,cG (1), Cauchy–Schwarz in (23) gives X h 2 i E[XJ2 ] ≤ Op,cH ,cG (1) E αη (cSjη )⊤ Mη (Pqη êJ ) . η∈Ω

Using (158) and the orthogonality of the exact-degree blocks, ℓ X

p X

∥cSj ∥22 = ∥c>p ∥22 ,

q=0

j=p+1

we conclude that

∥Pq êJ ∥22 = ∥êJ ∥22 ,

E[XJ2 ] ≤ Od n−1 ∥c>p ∥22 ∥êJ ∥22 = Od n−1 ∥c>p ∥22 .

Finally, the same bound holds for the terms containing ∆1 or ∆2 : these are simpler, since the extra factor ∆1 or ∆2 is controlled by its spectral norm, exactly as in the estimates (148). Then using hypercontractivity and Markov’s inequality we conclude the proof.

K.8

Proof of Proposition 33

Note that Sherman-Morrison-Woodbury formula for Kλ−1 gives −1 Kλ−1 = (ρ + λ)−1 H −1 − n1 (ρ + λ)−1 H −1 Φ≤p G−1 Φ⊤ . ≤p H

(159)

Then correspondingly, we have the decomposition −1 ⊤ 1 ⊤ n cSp+1 ΦSp+1 Kλ ΦSp+1 eJ 1 = (ρ+λ)

⊤ ⊤ −1 −1 −1 1 ⊤ ΦSp+1 eJ − n12 c⊤ Φ≤p G−1 Φ⊤ ΦSp+1 eJ Sp+1 ΦSp+1 H ≤p H n cSp+1 ΦSp+1 H

}

:=DJ

We now expand D as a whole. Since both summands in D involve the decomposition of H −1 , it is more convenient to classify terms only after expanding the entire expression. For brevity, write Φ := ΦSp+1 , c := cSp+1 . Using

H −1 = In + ∆H + ∆1 ,

G−1 = I + ∆G + ∆2 ,

we obtain DJ =

1 ⊤ ⊤ c Φ (In + ∆H + ∆1 )ΦeJ n 1 − 2 c⊤ Φ⊤ (In + ∆H + ∆1 )Φ≤p (I + ∆G + ∆2 )Φ⊤ ≤p (In + ∆H + ∆1 )ΦeJ . n

We decompose the full expansion of DJ into three disjoint classes: DJ = DJlead + DJcorr + DJrem . The leading part collects the terms with no ∆-factor: DJlead :=

1 ⊤ ⊤ 1 c Φ ΦeJ − c⊤ Φ⊤ Πp ΦeJ , n n

Πp := n1 Φ≤p Φ⊤ ≤p .

The structured correction part collects the terms that contain at least one factor ∆H or ∆G , but no factor ∆1 or ∆2 . Define Πp,G := n1 Φ≤p ∆G Φ⊤ ≤p 72

and

n o Mp := ∆H Πp , Πp ∆H , ∆H Πp ∆H , ∆H Πp,G , Πp,G ∆H , ∆H Πp,G ∆H , Πp,G .

Then we have the expression DJcorr :=

X 1 1 ⊤ ⊤ c Φ ∆H Φej − c⊤ Φ⊤ M ΦeJ . n n M ∈Mp

Finally, we define

DJrem := DJ − DJlead − DJcorr .

By construction, every term in DJrem contains at least one factor ∆1 or ∆2 . Thus, the three classes are: 1. DJrem : terms containing ∆1 or ∆2 ; 2. DJcorr : terms containing ∆H or ∆G , but not ∆1 or ∆2 ; 3. DJlead : terms containing no ∆-factor. We treat these three classes separately. An analogous analysis to equation (148) immediately gives

Analysis of DJrem .

|DJrem |2 2 √1 Φc n 2

≤

√1 Φc n

2 2 √1 ΦeJ ∥∆1 ∥op n 2 2 2 2 (2 H −1 op G−1 op + 2 2

= Od,P (d−p−δ log(d)) ∥c∥2

2 √1 Φ≤p n op

H −1 op )

2 √1 ΦeJ n 2

2 √1 ΦeJ . n 2

(160)

where we use the estimate of H −1 op = Od,P (1) (Eq. (141)), G−1 op = Od,P (1) (Eq. (143)),

√1 Φ≤p = n op

Od,P (1) (Lemma 7) and √1n Φc = Od,P (log(d)) ∥c∥2 (Proposition 26). 2 Then we take the expectation on both sides and obtain 2 2 E |DJrem |2 = Od,P (d−p−δ log2 (d)) ∥c∥2 ∥eJ ∥ . Analysis of DJlead .

(161)

Applying Lemma 39 immediately gives " 2 # 1 ⊤ ⊤ 2 2 E c Φ Φe)J = (c⊤ eJ )2 + Od (n−1 ) ∥c∥2 ∥eJ ∥2 , n

and " E

1 ⊤ ⊤ c Φ Πp ΦeJ n

2 #

= Od (d−2δ )(c⊤ eJ )2 + Od (n−1 ) ∥c∥2 ∥eJ ∥2 .

Therefore, we can deduce E

(DJlead − c⊤ eJ )2

" ≤ 2E

1 ⊤ ⊤ a Φ ΦeJ − (c⊤ eJ ) n

2 #

" + 2E 2

= Od (d−2δ )(c⊤ eJ )2 + Od (n−1 ) ∥c∥2 ∥eJ ∥2 . 73

1 ⊤ ⊤ c Φ Πp ΦeJ n

2 # (162)

Analysis of DJcorr . We use the following proposition to bound each summand in DJcorr . The proof appears in Appendix K.13. Proposition 37. Define Πp := n1 Φ≤p Φ⊤ ≤p ,

Πp,G := n1 Φ≤p ∆G Φ⊤ ≤p ,

and the finite collection n o Mp := ∆H Πp , Πp ∆H , ∆H Πp ∆H , ∆H Πp,G , Πp,G ∆H , ∆H Πp,G ∆H , Πp,G . Then, for every M ∈ Mp , the following holds " 2 # 1 ⊤ ⊤ 2 2 c Φ M Φ eJ = Od (d−2δ )(c⊤ eJ )2 + Od (n−1 ) ∥c∥2 ∥eJ ∥2 . E n2 Using Cauchy Schwartz inequality we obtain " 2 # X 1 ⊤ ⊤ 2 2 corr 2 (DJ ) ≤ |Mp | E c Φ M Φ eJ = Od (d−2δ )(c⊤ eJ )2 + Od (n−1 ) ∥c∥2 ∥eJ ∥2 . n2

(163)

(164)

M ∈Mp

Lastly, combining equation (162), (161) and (164) and using Cauchy-Schwartz inequality again gives E (D − c⊤ eJ )2 ≤ 3E (DJlead − c⊤ eJ )2 + 3E (DJrem )2 + 3E (DJcorr )2 2

= Od (d−2δ )(c⊤ eJ )2 + Od,P (d−p−δ ) log(d) ∥c∥2 ∥eJ ∥2 .

(165)

Using hypercontractivity, we have the following holds uniformly for all {u[j] }dj=1 2 2 E (D − c⊤ eJ )2 = Od (d−2δ+ϵ )(c⊤ eJ )2 + Od,P (d−p−δ+ϵ ) ∥c∥2 ∥eJ ∥2 .

(166)

With ∥eJ ∥2 = Od (d−p−1 ), taking the summation over J ∈ Sp and using Markov’s inequality completes the proof.

K.9

Proof of Proposition 34

Using the decomposition of Kλ−1 in (159), we obtain 1 ⊤ ⊤ −1 c Φ K ΦSp+1 eJ n Q Q λ 1 = c⊤ Φ⊤ H −1 ΦSp+1 eJ (ρ + λ)n Q Q 1 −1 − c⊤ Φ⊤ H −1 Φ≤p G−1 Φ⊤ ΦSp+1 eJ . ≤p H (ρ + λ)n2 Q Q

(167)

It therefore suffices to bound the two terms on the right-hand side. We expand H −1 and G−1 using (144) and (145), and decompose Φ≤p into exact-degree blocks exactly as in the proof of Proposition 33. In particular, every occurrence of an off-diagonal factor is expanded as off(B) = B − Diag(B). Thus each of the two terms in (167) becomes a finite sum of expressions of the form cη c⊤ Q Mη e J ,

where cη is deterministic, the number of summands is Op,cH ,cG (1), and Mη is a matrix chain of the form covered by Lemma 39. The left endpoint family of Mη is always Q, while the right endpoint family is always Sp+1 . Now the key simplification is that Q ∩ Sp+1 = ∅. Hence the endpoint families in every such chain are distinct, so the leading term in Lemma 39 vanishes automatically. Therefore, for every η, h 2 i E c⊤ = Od (nwη(1) wη(mη +1) )2 n−1 ∥cQ ∥22 ∥eJ ∥22 . Q Mη e J The deterministic coefficients cη are estimated exactly as in the proof of Proposition 33; in particular, |cη |2 (nwη(1) wη(mη +1) )2 n−1 = Od (n−1 ). Since there are only Op,cH ,cG (1) summands, we obtain " 2 # 1 ⊤ ⊤ −1 E cQ ΦQ H ΦSp+1 eJ = Od,P d−p−δ+ϵ ∥cQ ∥22 ∥eJ ∥22 , n

(168)

and similarly " E

1 ⊤ ⊤ −1 −1 c Φ H Φ≤p G−1 Φ⊤ ΦSp+1 eJ ≤p H n2 Q Q

2 #

(169)

= Od,P d−p−δ+ϵ ∥cQ ∥22 ∥eJ ∥22 .

Substituting (168) and (169) into (167) and using Cauchy-Schwartz inequality gives " 2 # 1 ⊤ ⊤ −1 cQ ΦQ Kλ ΦSp+1 eJ E = Od,P d−p−δ+ϵ ∥cQ ∥22 ∥eJ ∥22 . n

(170)

With ∥eJ ∥2 = Od (d−p−1 ), then taking the summation over J ∈ Sp and applying Markov’s inequality completes the proof.

K.10

Proof of Proposition 35

Note that for any ∆ = od,P (1), we can express (I + ∆)−1 using the Neumann series: (I + ∆)−1 = I +

∞ X

(−∆)k .

k=1

We compute the inverse of H (114) using Neumann series yielding q  2p+1 ∞ X X  −1  H −1 = In + − (ρ + λ)−1 dj θj d−j off(ΦSj Φ⊤ ∆ .  Sj ) − (ρ + λ) | {z } q=1 j=p+1

(171)

:=cj

Note that the triangle inequality on Eq. (31) and (32) shows for j ≥ p + 1 we have (172)

(p+δ−j)/2+ϵ d−j off(ΦSj Φ⊤ ) = od,P (1), Sj ) = Od,P (d

where off(·) denotes the off-diagonal component of a matrix. Next, we show that we can truncate 1

higher-order terms in H −1 at constant order to achieve a desired error on the order of Od,P (n− 2 ).

We set a constant cH := ⌈ p−1+2δ+2ϵ 1−δ−2ϵ ⌉. Using the triangle inequality and the formula for the summation of a geometric series, we can derive  q 2p+1 ∞ ∞ X X X −1   −cj d−j off(ΦSj Φ⊤ ) + (ρ + λ) ∆ ≤ (Od,P (d(δ−1)/2+ϵ ))q Sj q=cH +1

j=p+1

q=cH +1

= Od,P (d((δ−1)/2+ϵ)(cH +1) ) 1

= Od,P (n− 2 ). We expand the binomial for each q and bound the terms that contain ∆, which yields  q 2p+1 cH X X −1   −cj d−j off(ΦSj Φ⊤ ∆ H −1 = In + Sj ) + (ρ + λ) q=1

= In +

cH X

j=p+1



2p+1 X

 q=1

q (173)

 + ∆1 , −cj d−j off(ΦSj Φ⊤ Sj )

j=p+1 1

for some ∆1 satisfying ∥∆1 ∥op = Od,P (n− 2 ). We denote ∆H :=

cH X



2p+1 X

 q=1

q  , −cj d−j off(ΦSj Φ⊤ Sj )

j=p+1

1−δ

which satisfies ∥∆H ∥op = Od,P (d− 2 +ϵ ) with a quick calculation using Eq. (172). We now plug the expression (173) for H −1 into G in Eq. (118). We further substitute ⊤ Φ≤p Φ≤p off + IL(p) and we obtain n

G = (ρ + λ)Λ + I + off

Φ⊤ ≤p Φ≤p n

for

Φ⊤ Φ≤p ≤p + √ ∆H √ + ∆ 2 , n n

1 √ 2 with ∥∆2 ∥op ≤ ∥Φ≤p / n∥op ∥∆1 ∥op = Od,P (n− 2 ). Expressing G−1 using Neumann series again gives us

−1

=I+

∞ X

−(ρ + λ)Λ − off

r=1

Φ⊤ ≤p Φ≤p n

Φ⊤ Φ≤p ≤p − √ ∆H √ + ∆ 2 n n

!r ,

Similarly, we can truncate the infinite summation while maintaining the same error order. Note that the spectral norm of each item in the parentheses is of the order od,P (1), which can be deduced using equation (172) and the following estimates: ! δ Φ⊤ ≤p Φ≤p −δ ∥Λ∥op = Od (d ), off = Od,P (d− 2 +ϵ ), n op

Φ⊤ ≤p √

= Od,P (1),

∥∆2 ∥op = Od,P (n−1/2 ),

where particularly the second estimate can be deduced from Lemma 9 since diag 76

⊤

Φ≤p Φ≤p n

= In .

p+ϵ ⌉) Therefore, an analogous argument for H −1 shows that we can find a constant cG = max(cH , ⌈ δ−2ϵ 1

and a matrix ∆3 that satisfies ∥∆3 ∥op = Od,P (n− 2 ) such that the following holds G−1 = I +

cG X

−(ρ + λ)Λ−1 − off

r=1

Φ⊤ ≤p Φ≤p n

Φ⊤ Φ≤p ≤p − √ ∆H √ n n

!r (174)

+ ∆3 .

Then we conclude the proof for Eq. (145).

K.11

Proof of Proposition 36

To simplify the notation, write W := n−1/2 Φ≤p , Then we have

Vj := |Sj |−1/2 ΦSj

(p + 1 ≤ j ≤ 2p + 1).

1 ⊤ Φ Φ≤p = W ⊤ W = I + off(W ⊤ W ), n ≤p

and for p + 1 ≤ j ≤ 2p + 1, cj ′ ⊤ − j off(ΦSj Φ⊤ Sj ) = cj off(Vj Vj ), d

c′j := −cj

|Sj | = Od (1). dj

Hence every summand in ∆H is a finite product of matrices off(Vj Vj⊤ ) with deterministic coefficient of the order Od (1). We also write RG := AG − PG (Λ), and denote

B1 := off(W ⊤ W ),

B0 := −(ρ + λ)Λ, Then we have

∆G =

cG X

B2 := W ⊤ ∆H W.

(B0 + B1 + B2 )r ,

r=1

and consequently AH = B2 ,

AG = (I + B1 )∆G ,

AHG = B2 ∆G .

We now decompose the low-degree space according to exact degree. For 0 ≤ q ≤ p, let Wq := n−1/2 ΦSq , and let Pq : RL(p) → RL(p) be the orthogonal projector onto the coordinates indexed by Sq . Then W = [W0 , . . . , Wp ],

Wq = W Pq ,

p X

Pq ,

Λ=

q=0

p X

β q Pq ,

q=0

βq :=

dq . n

We have the decomposition B1 = off(W ⊤ W ) =

p X

off(Wq⊤ Wq ) +

0≤q,q ≤p q̸=q ′

Next, write X

cω Hω ,

Wq⊤ Wq′ .

′

q=0

∆H =

Hω :=

ω∈ΩH

mω Y ℓ=1

off Vjω,ℓ Vj⊤ , ω,ℓ

(175)

where ΩH is a finite index set and each cω = Od (1). Then we have B2 =

p X

(176)

cω Wq⊤ Hω Wq′ .

ω∈ΩH q,q ′ =0

Thus every word contributing to AH , RG , or AHG can be expanded using (175), (176) and B0 = −(ρ + λ)

p X

(177)

β q Pq .

q=0

Since p, cH , cG are fixed, this produces only Op,cH ,cG (1) summands. We also define the blockwise quantities Bt (b1 , b2 ) :=

p X

2 βq2t (Pq b1 )⊤ (Pq b2 ) ,

(178)

t ≥ 0.

q=0

Next we present the key claim. The proof appears in Appendix K.14. Proposition 38. Let T be one expanded word contributing to AH , RG , or AHG . Assume that T contains exactly s copies of B0 , and that T contains at least one off-diagonal factor, namely either a factor off(Wq⊤ Wq ) or one of the factors inside some Hω . Then s X 2 E (b⊤ T b ) ≤ Od (d−2δ ) Bt (b1 , b2 ) + Od n−1 ∥b1 ∥22 ∥b2 ∥22 . 2 1

(179)

t=0

We now return to AH , RG , and AHG . The pure Λ-words occur only when we take the identity part of W ⊤ W and every factor in ∆G is B0 ; collecting these terms gives precisely PG (Λ) =

cG X

r −(ρ + λ)Λ .

r=1

Hence, by definition, every word contributing to RG contains at least one off-diagonal factor. The same is clearly true for every word in AH = B2 and in AHG = B2 ∆G . Applying the claim word-by-word, and using that each word contains at most cG copies of B0 , we obtain, for each X ∈ {AH , RG , AHG }, cG X 2 Od (d−2δ ) Bt (b1 , b2 ) + Od n−1 ∥b1 ∥22 ∥b2 ∥22 . E (b⊤ ≤ 1 Xb2 )

(180)

t=0

The final piece of the proof is the upper bound for Eq. (180) completes the proof.

t=0 Bt (b1 , b2 ). Combining the bound below with

Claim 5. The following estimate holds: cG X

! 2

Bt (c≤p , êJ ) = Od (1) (c≤p êJ ) + Od (d

−4

)·

t=0

K.12

c2J\{i}

(181)

i∈J

Proof of Claim 4

Expanding PG (Λ) gives ⊤ c⊤ ≤p PG (Λ)êJ = c≤p

cG X

! −(ρ + λ)Λ

êJ

r=1

cG X

r (−(ρ + λ))r c⊤ ≤p Λ êJ .

r=1

(182)

Talking the square and applying Cauchy–Schwarz inequality yields c⊤ ≤p PG (Λ)êJ

≤ cG

cG X

r (ρ + λ)2r c⊤ ≤p Λ êJ

r=1

For each r = 1, . . . , cG , we have r c⊤ ≤p Λ êJ

2 −r −1 = c⊤ D E:L(p),J ≤p (nD)  2 X X −1 ′ = (nDJ⊔{i} )−r cJ⊔{i} ui + (nDJ\{i} )−r d−1 DJ,J DJ\{i} cJ\{i} ui  i∈S

i∈[d]\J

2

 ≤ Od (d−2r(p+δ−|J|−1) ) 

cJ⊔{i} ui 

i∈[d]\J

}

:=I 2

!2 −2r(p+δ−|J|+1)

+ Od (d

) d

−1

−1 ′ DJ,J DJ\{i}

cJ\{i} ui

i∈J

}

:=J 2

where we use the basic inequality (a + b)2 ≤ 2(a2 + b2 ) in the last inequality. Now using Young’s inequality ax2 + by 2 ≤ 2a(x + y)2 + (2a + b)y 2 , we have Od (d−2r(p+δ−|J|−1) )I 2 + Od (d−2r(p+δ−|J|+1) )J 2 ≤ Od (d−2r(p+δ−|J|−1) )(I + J )2 + Od (d−2r(p+δ−|J|+1) )J 2 ! 2

= Od (d−2r(p+δ−|J|−1) ) (c≤p êJ ) + Od (d−2r(p+δ−|J|+1) )Od (d−4 ) ·

c2J\{i}

i∈J

where we apply the bound for J 2 from equation (195) in the last equality. When r = 1 and |J| = p − 1, we have the smallest upper bound which is ! X 2 −2δ −8−2δ 2 Od (d ) (c≤p êJ ) + Od (d ) cJ\{i} .

(183)

i∈J

Then we complete the proof.

K.13

Proof of Proposition 37

Set V := Vp+1 = |Sp+1 |−1/2 ΦSp+1 , and define QM :=

1 ⊤ ⊤ a ΦSp+1 M ΦSp+1 e, n

γ :=

|Sp+1 | , n

M ∈ Np .

Since ΦSp+1 = |Sp+1 |1/2 V,

Πp =

1 ⊤ Φ≤p Φ⊤ ≤p = W W , n

Πp,G =

1 ⊤ Φ≤p ∆G Φ⊤ ≤p = W ∆G W , n

every expanded word contributing to QM has the form γ a⊤ V ⊤ H0 W Ξ W ⊤ H1 V e, 79

(184)

where H0 , H1 are either I or one of the words appearing in the expansion of ∆H , and Ξ is either I or one of the words appearing in the expansion of ∆G . As in the low-degree endpoint proof, we decompose the low-degree space by exact degree: for 0 ≤ q ≤ p, Wq := W Pq = n−1/2 ΦSq ,

Λ=

p X

βq P q ,

dq . n

βq :=

q=0

Expanding every occurrence of B0 , B1 , B2 , and then resolving every occurrence of W and W ⊤ into the exact-degree blocks Wq , Wq⊤ , we obtain a finite decomposition QM =

(185)

γ αη a⊤ Mη e,

η∈Ω

where |Ω| = Op,cH ,cG (1), and for each η ∈ Ω: 1. αη is the product of all deterministic scalar coefficients coming from the chosen words in ∆H and ∆G , together with the factors βq contributed by the selected copies of B0 . In particular, |αη | = Od (1). 2. Mη is the corresponding normalized matrix chain obtained after removing those scalar coefficients. Concretely, Mη is a product of matrices chosen from V ⊤,

Wq⊤ ,

Wq ,

off(Vj Vj⊤ ),

off(Wq⊤ Wq ),

Wq⊤ Wq′ ,

with endpoint family Sp+1 on both sides. We claim that every summand in (185) satisfies h 2 i E γ αη a⊤ Mη e ≤ od (1) (a⊤ e)2 + Od n−1 ∥a∥22 ∥e∥22 .

(186)

There are two cases. Case 1: Mη contains at least one off-diagonal replacement. Let Sη1 , . . . , Sηmη +1 be the sequence of families appearing in the chain Mη , and let Vη , Eη be the corresponding sets of samplespace and feature-space off-diagonal replacements. Then Mη satisfies the hypotheses of Corollary 42. The crucial point is that Mη always contains at least one internal low-degree family Sq : indeed, every word in (184) contains the middle factor W Ξ W ⊤ , and after exact-degree expansion this contributes at least one block Wq or Wq⊤ . Therefore κMη :=

|Sηk | |Sq | ≤ = Od (d−2δ ), 2≤k≤mη n n min

since p is fixed and |Sq | ≍ dq = Od (d−δ ) for all q ≤ p. Applying Corollary 42 to Mη , we obtain E (a⊤ Mη e)2 ≤ Od κMη |Eη | (nw(1) w(mη +1) )2 (a⊤ e)2 + Od (nw(1) w(mη +1) )2 n−1 ∥a∥22 ∥e∥22 . Here the endpoint family is Sp+1 on both sides, so w

(1)

(mη +1)

−1/2

= |Sp+1 |

(nw 80

(1)

(mη +1) 2

) =

n |Sp+1 |

= γ −2 .

(187)

Also |Eη | = Op,cH ,cG (1), since the total number of factors is uniformly bounded. Hence (187) simplifies to E (a⊤ Mη e)2 ≤ Od (d−2δ ) γ −2 (a⊤ e)2 + Od γ −2 n−1 ∥a∥22 ∥e∥22 . Multiplying by γ 2 |αη |2 , and using |αη | = Od (1), yields (186). Case 2: Mη contains no off-diagonal replacement. This can happen only when M = Πp,G , H0 = H1 = I, and the chosen word in ∆G is a pure B0 -word. Indeed, every occurrence of ∆H contributes an off-diagonal factor, and every occurrence of B1 does as well. Thus for some 1 ≤ r ≤ cG ,

Ξ = B0r and (185) reduces to QM = γ

p X

αq a⊤ V ⊤ Wq Wq⊤ V e,

αq = Od (βqr ).

(188)

q=0

For each fixed q, the chain V ⊤ Wq Wq⊤ V has endpoint family Sp+1 on both sides and only one low-degree family Dq . Therefore Lemma 39 gives h 2 i E a⊤ V ⊤ Wq Wq⊤ V e ≤ Cd γ −2 (a⊤ e)2 + n−1 ∥a∥22 ∥e∥22 . Using Cauchy–Schwarz in (188), we obtain E[Q2M ] ≤ Cp

p X

h 2 i |αq |2 γ 2 E a⊤ V ⊤ Wq Wq⊤ V e .

q=0

Hence we have

p X

E[Q2M ] ≤ Cd

! |αq |

(a⊤ e)2 + n−1 ∥a∥22 ∥e∥22 .

q=0

Since |αq | = Od (βqr ), r ≥ 1, and max βq =

0≤q≤p

we have

p X

dp = Od (d−δ ) n

|αq |2 = Od (d−2δ ),

q=0

so (186) follows in this case as well. Combining the two cases proves (186). Finally, since |Ω| = Op,cH ,cG (1), Cauchy–Schwarz in (185) yields, for every M ∈ Np , " 2 # 1 ⊤ ⊤ E a ΦSp+1 M ΦSp+1 e ≤ Od (d−2δ ) (a⊤ e)2 + Od n−1 ∥a∥22 ∥e∥22 . n This proves Proposition 37.

K.14

Proof of Proposition 38

Expanding T using (175), (176), and the block decomposition of B0 (177), we may write X ⊤ η b ) M (P η b ) , b⊤ αη (Pq− 1 η q+ 2 1 T b2 = η∈Ω(T )

(189)

where Ω(T ) is a finite set of cardinality Op,cH ,cG (1), each coefficient αη is deterministic and has the form s Y αη = Od βuη,ℓ , uη,ℓ ∈ {0, . . . , p}, (190) ℓ=1

and each Mη is a normalized matrix chain whose families belong to S0 , . . . , S2p+1 with the same off-diagonal replacements as in the original word T . In particular, Lemma 39 and Corollary 42 apply to each Mη . Fix one η ∈ Ω(T ). η η If q− ̸= q+ , then the endpoint families differ, and Lemma 39 gives h 2 i ⊤ 2 2 η b ) M (P η b ) η b ∥ ∥P η b ∥ E (Pq− = Od n−1 ∥Pq− (191) 1 η 1 2 q+ 2 q+ 2 2 . η η Now suppose q− = q+ = q. If every low-degree family occurring in Mη is equal to Sq , then Mη has the same low-degree family at both endpoints and, by assumption, contains at least one off-diagonal replacement. Hence Corollary 42 applies and yields h 2 i 2 E (Pq b1 )⊤ Mη (Pq b2 ) = Od (d2(q−p−δ) ) (Pq b1 )⊤ (Pq b2 ) + Od n−1 ∥Pq b1 ∥22 ∥Pq b2 ∥22 .

If Mη contains another low-degree family Sq′ with q ′ ̸= q, then Mη contains at least two distinct families of size o(n), namely Sq and Sq′ . Therefore we are in the second case of Lemma 39, and h 2 i 2 ′ E (Pq b1 )⊤ Mη (Pq b2 ) = Od (d2(q ∧q)−2p−2δ ) (Pq b1 )⊤ (Pq b2 ) + Od n−1 ∥Pq b1 ∥22 ∥Pq b2 ∥22 . Thus, whenever the endpoints agree, we have the uniform bound h 2 i 2 E (Pq b1 )⊤ Mη (Pq b2 ) ≤ Od (d2(q−p−δ) ) (Pq b1 )⊤ (Pq b2 ) + Od n−1 ∥Pq b1 ∥22 ∥Pq b2 ∥22 .

(192)

It remains to account for the coefficient αη . If the endpoints agree and are equal to q, let tη (q) ∈ {0, . . . , s} denote the number of copies of B0 that act on the degree-q block in the summand indexed by η. Then, by Eq. (190), we have |αη |2 ≲d βq2tη (q) βp2(s−tη (q)) .

(193)

2(s−t (q))

η Since βp = dp /n = od (1), the factor βp is od (1) whenever tη (q) < s, while if tη (q) = s it is simply 1. Combining (192) with (193), and using (191) when the endpoints differ, we obtain

p s h X 2 i X 2 ⊤ η b ) M (P η b ) E αη (Pq− ≤ Od (d−2δ ) βq2t (Pq b1 )⊤ (Pq b2 ) + Od n−1 ∥b1 ∥22 ∥b2 ∥22 . 1 η q+ 2 t=0

q=0

This is exactly s h 2 i X ⊤ η b ) M (P η b ) E αη (Pq− ≤ Od (d−2δ ) Bt (b1 , b2 ) + Od n−1 ∥b1 ∥22 ∥b2 ∥22 . 1 η q+ 2 t=0

Since |Ω(T )| = Op,cH ,cG (1), Cauchy–Schwarz in (189) proves (179).

K.15

Proof of Claim 5

q−p−δ Since (nD)−1 ) = od (1) for q ≤ p, expanding Sq is a diagonal matrix with each entry of the order Od (d Bt gives p D cG cG X E2 X X Bt (c≤p , êJ ) = cSq , (nDSq )−2t DS−1 ESq ,J q t=0

t=0 q=0

= (1 + od (1))

p D E2 X cSq , DS−1 E . Sq ,J q q=0

(194)

Plugging the expression of DS−1 ESq ,J namely (108) into the above equation gives q  2 !2 p D E2 X X X −1 −1 −1 ′ cSq , DSq ESq ,J =  cJ\{i} ui . cJ⊔{i} ui  + d DJ,J DJ\{i} q=0

i∈J

i∈[d]\J

:=I 2

}

:=J 2

By Cauchy-Schwarz, we obtain ! −1 ′ J 2 ≤ d−2 · (DJ,J DJ\{i} )2 · (

u2i )

i∈J

(195)

c2J\{i}

i∈J

! = Od (d

−4

)·

c2J\{i}

(196)

i∈J

thus by Young’s inequality we have ! 2

I + J ≤ 2(I + J ) + 3J = 2 (c≤p êJ ) + Od (d

−4

)·

c2J\{i}

i∈J

Combining it with Eq. (194) completes the proof.

Product of random Fourier-Walsh matrices

In this section, we state and prove a key lemma for our main results. The proof is adapted from [57]. Lemma 39. Fix families S1 , . . . , Sm+1 ⊆ 2[d] , and let ⊤ ⊤ M := A⊤ 1 A2 A2 · · · Am Am Am+1 ,

Ai := w(i) XSi ,

where w(i) = n−1/2 ∧ |Si |−1/2 . Assume the following: (i) every S ∈ Si satisfies |S| ≤ pi ; (ii) if Si ∩ Sj ̸= ∅, then Si = Sj . Let b1 ∈ R|S1 | , b2 ∈ R|Sm+1 | , and set

κ := (n w(1) w(m+1) )2 .

If S1 ̸= Sm+1 , then the following estimate for second-moment holds: h 2 i E b⊤ = κ · Od n−1 ∥b1 ∥22 ∥b2 ∥22 . 1 M b2 Assume now that S1 = Sm+1 . If every distinct internal family Sk ̸= S1 , 2 ≤ k ≤ m, satisfies n = O(|Sk |), then the second-moment expansion h 2 i 2 −1 E b⊤ = κ (b⊤ ∥b1 ∥22 ∥b2 ∥22 (197) 1 M b2 1 b2 ) + O d n holds, and the estimate for first-moment expansion also holds: 1 √ −2 ⊤ E b⊤ M b = κ b b + O n ∥b ∥ ∥b ∥ . 2 d 1 2 2 2 1 1 2

(198)

In the complementary case, which necessarily has m ≥ 2, the following estimates hold: h 2 i |Sk |2 ⊤ 2 −1 2 2 ⊤ (b1 b2 ) + Od n ∥b1 ∥2 ∥b2 ∥2 E b1 M b 2 = κ Od min 2≤k≤m n2

(199)

and 1 ⊤ √ |Sk | −2 ⊤ E b1 M b2 = κ Od min (b1 b2 ) + Od n ∥b1 ∥2 ∥b2 ∥2 . 2≤k≤m n

(200)

Proof of Lemma 39. We first prove the second-moment estimates. The first-moment estimates follow from the same collapse procedure with one copy of the chain instead of two; the necessary modifications are given at the end of the proof. For readability, introduce the random variable ⊤ ⊤ ⊤ S := b⊤ 1 A1 (A2 A2 ) · · · (Am Am )Am+1 b2

Expanding the square gives the identity " X ⊤ E[S] = E (A1 b1 )i1 (A2 A⊤ 2 )i1 ,i2 · · · (Am Am )im−1 ,im (Am+1 b2 )im i1 ,...,i2m ∈[n]

# .

⊤ × (A1 b1 )im+1 (A2 A⊤ 2 )im+1 ,im+2 · · · (Am Am )i2m−1 ,i2m (Am+1 b2 )i2m

For each partition π ∈ Π2m , define the exact equality class Iπ := {(i1 , . . . , i2m ) ∈ [n]2m : is = it if and only if s ∼π t}. The resulting partition decomposition is E[S] =

E[Sπ ],

π∈Π2m

where the contribution associated with π is defined by " X ⊤ E (A1 b1 )i1 (A2 A⊤ E[Sπ ] := 2 )i1 ,i2 · · · (Am Am )im−1 ,im (Am+1 b2 )im (i1 ,...,i2m )∈Iπ

# .

⊤ × (A1 b1 )im+1 (A2 A⊤ 2 )im+1 ,im+2 · · · (Am Am )i2m−1 ,i2m (Am+1 b2 )i2m

We now fix π ∈ Π2m and estimate E[Sπ ]. At each stage of the reduction, the surviving sample positions are relabeled from left to right as i1 , . . . , iq . The current product consists of factors of the form (A1 b1 )ir ,

(Au A⊤ u )ir ,ir+1 ,

(Am+1 b2 )ir ,

with family labels inherited from the unreduced product. We repeatedly use the following two reductions. First, suppose that the current partition forces two adjacent sample positions to be equal, say ir = ir+1 , and that the factor between them is (Au A⊤ u )ir ,ir+1 . The diagonal identity gives the scalar value (u) 2 (Au A⊤ ) |Su | = 1 ∧ u )ir ,ir = (w

|Su | . n

Thus the Gram factor is replaced by a scalar, which we denote by CF (u) := (w(u) )2 |Su | = 1 ∧

|Su | . n

One of the two equal sample variables is then deleted. We call this operation a feature collapse. 84

Second, suppose that {r} is a singleton block of the current partition. For fixed values of the other current block labels, the label ir has n + Om (1) admissible choices. Define the endpoint pairing by X ⟨b1 , b2 ⟩∩ := (b1 )S (b2 )S . S∈S1 ∩Sm+1

Also set

CS (u, v) := (n + Om (1))w(u) w(v) 1Su =Sv .

Independence across samples and orthogonality of distinct feature families give the following local replacements: ⊤ ⊤ (Au A⊤ u )ir−1 ,ir (Av Av )ir ,ir+1 ⇝ CS (u, v)(Au Au )ir−1 ,ir+1 , (A1 b1 )ir (Av A⊤ v )ir ,ir+1 ⇝ CS (1, v)(A1 b1 )ir+1 , (Au A⊤ u )ir−1 ,ir (Am+1 b2 )ir ⇝ CS (u, m + 1)(Am+1 b2 )ir−1 , (A1 b1 )ir (Am+1 b2 )ir ⇝ CS (1, m + 1)⟨b1 , b2 ⟩∩ . Here ⇝ means that the contribution is unchanged after summing over the singleton sample label and deleting that variable. In the first line, the family label on the remaining Gram factor is immaterial because the factor 1Su =Sv forces Su = Sv . We call these operations sample collapses. We iterate the two reductions until neither applies. Let π ♯ be the final partition, and let cπ be the product of all scalar coefficients generated during the reduction. With the convention that any endpoint pairing ⟨b1 , b2 ⟩∩ produced before the procedure stops is kept as an external scalar factor, the collapse procedure gives the relation  cπ ⟨b1 , b2 ⟩2∩ , |π ♯ | = 0, E[Sπ ] = (201) c E[S ♯ ], |π ♯ | ≥ 1. π π All feature-collapse and sample-collapse coefficients have absolute value at most one, so the coefficient cπ satisfies 0 ≤ cπ ≤ 1. We now record the coefficient estimate for fully collapsed partitions. The proof appears in Appendix L.1. Proposition 40 (Fully collapsed coefficient). Assume S1 = Sm+1 . Fix π ∈ Π2m such that |π ♯ | = 0 and cπ > 0. For a distinct family T , write wT := n−1/2 ∧ |T |−1/2 . Let NF (T ) be the number of feature collapses performed on T , and let NS (T ) be the number of nonterminal sample collapses performed on T , excluding the two terminal endpoint collapses that produce the two factors ⟨b1 , b2 ⟩∩ . Then the coefficient admits the factorization N (T ) NS (T ) Y |T | F n −1 −1 κ cπ = 1 + Om (n ) 1∧ 1∧ . n |T | T

In particular, if no collapse produces a loss, then the coefficient satisfies cπ = κ 1 + Om (n−1 ) . Moreover, if at least one collapse on a family T produces a loss, then |T | n cπ = κ · O min , . n |T | Finally, if T ̸= S1 is a distinct internal family, say T = Sk for some 2 ≤ k ≤ m, and |T | < n, then every fully collapsed second-moment contribution satisfies 2 |T | cπ = κ · O . n2 85

We next bound the residual terms, namely the terms for which the collapse procedure stops before the entire product disappears. The proof appears in Appendix L.2. Proposition 41. Fix π ∈ Π2m , and suppose that the recursive collapse procedure stops at a final partition π ♯ with |π ♯ | ≥ 1. Then the following estimate for residual holds: E[Sπ♯ ] = Od (n−1 )∥b1 ∥22 ∥b2 ∥22 We now assemble the second-moment estimates. The bound κ−1 = Od (1) follows from |S1 |, |Sm+1 | ≤ 2d . Proposition 41, together with 0 ≤ cπ ≤ 1, implies that every non-fully-collapsed partition contributes a term satisfying cπ E[Sπ♯ ] = κ · Od n−1 ∥b1 ∥22 ∥b2 ∥22 . (202) If S1 ̸= Sm+1 , then the endpoint pairing vanishes: ⟨b1 , b2 ⟩∩ = 0. Thus every fully collapsed contribution vanishes. Since the number of partitions of [2m] depends only on m, the non-fully-collapsed bound in (202) gives h 2 i E b⊤ M b = κ · Od n−1 ∥b1 ∥22 ∥b2 ∥22 . 2 1 Assume now that S1 = Sm+1 holds. The endpoint pairing is then ⟨b1 , b2 ⟩∩ = b⊤ 1 b2 . Suppose first that every distinct internal family Sk = ̸ S1 , 2 ≤ k ≤ m, satisfies n = O(|Sk |). The fully collapsed contribution in which all collapses are loss-free has coefficient κ 1 + Om (n−1 ) . Every other fully collapsed contribution contains at least one lossy collapse. If the loss occurs on a family T with |T | < n, then the loss |T |/n is Od (n−1 ). If the loss occurs on a distinct internal family T ̸= S1 with |T | > n, then the loss is n/|T |; the assumption n = O(|T |), together with |T | ≤ 2d , implies that this loss is also Od (n−1 ). Hence all fully collapsed contributions except the loss-free one are absorbed into the Od (n−1 ) error. The fully collapsed terms therefore satisfy X 2 −1 E[Sπ ] = κ (b⊤ ∥b1 ∥22 ∥b2 ∥22 . 1 b2 ) + Od n π∈Π2m |π ♯ |=0

Combining this estimate with (202) gives h 2 i 2 −1 E b⊤ = κ (b⊤ ∥b1 ∥22 ∥b2 ∥22 , 1 M b2 1 b2 ) + O d n which proves (197). It remains to consider the complementary regime. In this case, there is a distinct internal family T = Sk ̸= S1 , 2 ≤ k ≤ m, with |T | < n along the relevant range. Define ρ2 :=

min

2≤k≤m Sk ̸=S1 , |Sk |<n

|Sk |2 . n2

By Proposition 40, every fully collapsed second-moment contribution must incur two feature-collapse losses from such an internal family. This gives the bound X 2 E[Sπ ] = κ · Od (ρ2 )(b⊤ 1 b2 ) . π∈Π2m |π ♯ |=0

Since all families are nonempty and have size at most 2d , replacing ρ2 by min2≤k≤m |Sk |2 /n2 changes only the Od (·) constant. After adding the non-fully-collapsed contribution from (202), we obtain h 2 i |Sk |2 ⊤ 2 −1 2 2 E b⊤ M b = κ O min (b b ) + O n ∥b ∥ ∥b ∥ , 2 d d 1 2 2 2 1 1 2 2≤k≤m n2 which proves (199). We finally prove the first-moment estimates. Set T := b⊤ 1 M b2 . The partition decomposition is now taken over [m], because there is only one copy of the chain. The same feature and sample collapses apply. A fully collapsed first-moment term has one terminal endpoint collapse instead of two, so the terminal coefficient has the form √ (n + Om (1))w(1) w(m+1) = κ 1 + Om (n−1 ) . The analogue of Proposition 41, proved by the same incidence-graph argument with one chain instead of two, gives the stronger residual estimate E[Tπ♯ ] = Od (n−1 )∥b1 ∥2 ∥b2 ∥2 for every non-fully-collapsed first-moment residual. Since κ−1/2 = Od (1), the weaker estimate needed below also holds: √ E[Tπ♯ ] = κ Od (n−1/2 )∥b1 ∥2 ∥b2 ∥2 . (203) Assume first that every distinct internal family Sk = ̸ S1 , 2 ≤ k ≤ m, satisfies n = O(|Sk |). The fully collapsed loss-free contribution has value √ κ 1 + Om (n−1 ) b⊤ 1 b2 . As in the second-moment argument, all other fully collapsed contributions are absorbed into the Od (n−1 ) error. The Cauchy–Schwarz inequality gives the auxiliary bound −1/2 n−1 |b⊤ ∥b1 ∥2 ∥b2 ∥2 . 1 b2 | ≤ n

Combining the fully collapsed estimate with (203) gives i √ h ⊤ E b⊤ κ b1 b2 + Od n−1/2 ∥b1 ∥2 ∥b2 ∥2 , 1 M b2 = which proves (198). In the complementary regime, define ρ1 :=

min

2≤k≤m Sk ̸=S1 , |Sk |<n

|Sk | . n

A fully collapsed first-moment contribution must incur one feature-collapse loss from a distinct internal family with size smaller than n. This gives the fully collapsed bound X √ E[Tπ ] = κ Od (ρ1 )(b⊤ 1 b2 ). π∈Πm |π ♯ |=0

Again, because all family sizes lie between 1 and 2d , replacing ρ1 by min2≤k≤m |Sk |/n changes only the Od (·) constant. Adding the non-fully-collapsed first-moment contribution from (203) yields √ |Sk | ⊤ −1/2 E b⊤ M b = κ O min (b b ) + O n ∥b ∥ ∥b ∥ , 2 d d 1 2 2 2 1 1 2 2≤k≤m n which proves (200). 87

Corollary 42. Fix families S1 , . . . , Sm+1 ⊆ 2[d] , and define ⊤ ⊤ M := A⊤ 1 A2 A2 · · · Am Am Am+1 ,

Ai := w(i) XSi ,

where the weights are given by w(i) := n−1/2 ∧ |Si |−1/2 . Assume that the families S1 , . . . , Sm+1 satisfy the hypotheses of Lemma 39. Assume moreover that the endpoint families agree, namely, S1 = Sm+1 . Let b1 , b2 ∈ R|S1 | and let V, E ⊆ {2, . . . , m} be such that V ∪ E ̸= ∅, V ∩ E ∪ (E + 1) = ∅,

E + 1 := {j + 1 : j ∈ E}.

Assume also that the adjacent matrices agree along E, that is, Aj = Aj+1

for every j ∈ E.

for every j ∈ V,

|Sj | ≤ n

Finally, assume the size conditions n ≤ |Sj |

for every j ∈ E.

Define MV,E to be the matrix obtained from M by the following replacements. For each j ∈ V, replace Aj A⊤ j

off(Aj A⊤ j ).

A⊤ j Aj+1

off(A⊤ j Aj+1 ).

Also, for each j ∈ E, replace Set κ := (n w(1) w(m+1) )2 . Let B denote the set of distinct internal families T ̸= S1 for which the first-regime condition n = O(|T |) fails. Equivalently, B is the set of distinct internal families that force the “otherwise” case of Lemma 39. Define  B = ∅, 0, ρ := |T |  min , B ̸= ∅. T ∈B n Then

h 2 i 2 −1 E b⊤ = κ · Od (ρ2 )(b⊤ ∥b1 ∥22 ∥b2 ∥22 . 1 MV,E b2 1 b2 ) + κ · Od n

Proof of Corollary 42. We first record the diagonal identities used in the expansion. If j ∈ V, then n ≤ |Sj |, and therefore w(j) = |Sj |−1/2 . Hence (j) 2 Diag(Aj A⊤ ) |Sj | In = In . j ) = (w

Similarly, if j ∈ E, then Aj = Aj+1 and |Sj | ≤ n. Thus w(j) = w(j+1) = n−1/2 , and hence (j) (j+1) Diag(A⊤ w I|Sj | = I|Sj | . j Aj+1 ) = n w

Consequently, for j ∈ V, we have Likewise, for j ∈ E, we have

⊤ off(Aj A⊤ j ) = A j Aj − I n . ⊤ off(A⊤ j Aj+1 ) = Aj Aj+1 − I|Sj | .

We now expand the off-diagonal constraints by inclusion–exclusion. For J ⊆ V and K ⊆ E, define the remaining index set by ΓJ,K := [m + 1] \ J ∪ (K + 1) = {i1 < · · · < ir }. On this reduced index set, define the corresponding chain by ⊤ ⊤ MJ,K := (A⊤ i1 Ai2 )(Ai2 Ai3 ) · · · (Air−1 Air ).

The condition

V ∩ E ∪ (E + 1) = ∅

ensures that the deletions coming from J and K do not interfere with one another. Therefore, inclusion– exclusion gives X X MV,E = (−1)|J|+|K| MJ,K . (204) J⊆V K⊆E

We next verify that each reduced chain MJ,K is still covered by Lemma 39. Deleting an index from V removes only a family Sj with n ≤ |Sj |, and therefore it cannot remove any family in B. Deleting j + 1 with j ∈ E only shortens a block of equal adjacent matrices, since Aj = Aj+1 . Thus no distinct bad internal family is removed entirely. The endpoint families are also unchanged. The left endpoint is never deleted. If the right endpoint is deleted, then K contains a terminal interval {ℓ, ℓ + 1, . . . , m} for some ℓ. In that case, the new right endpoint is ℓ, and the equalities along this terminal block imply Aℓ = Aℓ+1 = · · · = Am+1 . Thus the right endpoint family and the right endpoint weight remain those of Sm+1 . Hence every reduced chain has the same endpoint families and the same value of κ. For notational convenience, write the index set of reduced chains as A := {(J, K) : J ⊆ V, K ⊆ E}. For a = (J, K) ∈ A, define

σa := (−1)|J|+|K| ,

Finally, denote

β := b⊤ 1 b2 ,

Za := b⊤ 1 MJ,K b2 . N := ∥b1 ∥2 ∥b2 ∥2 .

With this notation, equation (204) becomes b⊤ 1 MV,E b2 =

σa Za .

(205)

a∈A

We now distinguish the two regimes of Lemma 39. First regime: B = ∅. In this case, every reduced chain lies in the first endpoint-matched regime of Lemma 39. We use the following mixed second-moment estimate: for all a, a′ ∈ A, E[Za Za′ ] = κ β 2 + Od n−1 N 2 . (206) This is the polarized form of the second-moment estimate in Lemma 39. Equivalently, it follows by applying the same collapse argument to the product Za Za′ . The fully collapsed loss-free contribution gives the common leading term κβ 2 , while all lossy or non-fully-collapsed contributions are absorbed into the Od (n−1 N 2 ) error. Since V ∪ E ̸= ∅, the alternating coefficients satisfy    X X X σa =  (−1)|J|   (−1)|K|  = 0. (207) a∈A

J⊆V

K⊆E

Using (205), (206), and (207), we obtain h X 2 i E b⊤ = σa σa′ E[Za Za′ ] 1 MV,E b2 a,a′ ∈A

!2 = κβ

+ κ · Od n−1 N 2

σa

a∈A −1

= κ · Od n

N2 .

Since ρ = 0 in the present regime, this gives the desired bound. Second regime: B ̸= ∅. In this case, every reduced chain remains in the second regime of Lemma 39. Indeed, as observed above, the deletions do not remove any family in B. Therefore, Lemma 39 gives, for every a ∈ A, E[Za2 ] = κ Od (ρ2 )β 2 + Od n−1 N 2 . (208) Since |A| ≤ 2m , Cauchy–Schwarz applied to the finite sum in (205) gives h X 2 i E b⊤ M b ≤ |A| E[Za2 ]. V,E 2 1 a∈A

Combining this estimate with (208), we obtain h 2 i E b⊤ = κ · Od (ρ2 )β 2 + κ · Od n−1 N 2 . 1 MV,E b2 Finally, substituting back β = b⊤ 1 b2 and N = ∥b1 ∥2 ∥b2 ∥2 , yields h 2 i 2 −1 E b⊤ = κ · Od (ρ2 )(b⊤ ∥b1 ∥22 ∥b2 ∥22 . 1 MV,E b2 1 b2 ) + κ · Od n This completes the proof.

L.1

Proof of Proposition 40

A feature collapse on T has contribution wT2 |T | = 1 ∧

|T | . n

Similarly, the contribution of a nonterminal sample collapse on T is n (n + Om (1))wT2 = 1 + Om (n−1 ) 1 ∧ . |T | The two terminal endpoint collapses have combined contribution 2 (n + Om (1))w(1) w(m+1) = κ 1 + Om (n−1 ) . Multiplying all collapse contributions gives the factorization N (T ) NS (T ) Y n |T | F 1∧ . κ−1 cπ = 1 + Om (n−1 ) 1∧ n |T |

(209)

Separating the feature families according to whether |T | < n or |T | > n, and noting that families with |T | = n contribute a factor equal to one, gives the equivalent representation κ−1 cπ = 1 + Om (n−1 )

Y |T |<n

|T | n

NF (T ) Y |T |>n

n |T |

NS (T ) .

(210)

The first two consequences follow immediately from (210). It remains to prove the last claim. Let T ̸= S1 be an internal feature family with |T | < n. In each copy of the chain, sample collapses can merge adjacent T -factors, but they cannot remove the final surviving T -Gram factor, because the endpoints have family S1 ̸= T . Consequently, each copy must contain at least one feature collapse on T . Since there are two copies of the chain, the bound NF (T ) ≥ 2 holds. Applying this bound to (210) gives cπ = κ · O

|T |2 n2

This proves the proposition.

L.2

Proof of Proposition 41

Let q be the number of surviving sample positions after the collapse procedure stops. Set r := |π ♯ |, so that r is the number of blocks in the final partition. Since no further sample collapse is possible, the partition π ♯ has no singleton block. This observation gives the counting inequality (211)

q ≥ 2r.

Moreover, since |π ♯ | ≥ 1, the integer r satisfies r ≥ 1. During the previous reductions, one entire copy of the chain may already have collapsed and produced a scalar endpoint factor ⟨b1 , b2 ⟩∩ . Let Θ denote the product of any such external endpoint factors, with the convention that Θ = 1 if no such factor is present. We keep Θ outside the graph notation below. Let G♯ be the incidence multigraph whose vertices are the blocks of π ♯ . Each surviving factor (Au A⊤ u )is ,it gives an internal edge e joining the two blocks containing s and t. Each surviving endpoint factor (A1 b1 )is or (Am+1 b2 )is gives a half-edge h attached to the block containing s. Since no further feature collapse is possible, no internal edge is a loop. Write E = E(G♯ ) and H = H(G♯ ). For each internal edge e ∈ E, let Fe = Sue be the corresponding feature family. The edge factor associated with e admits the representation X (Aue A⊤ xSi e xSj e , ωe := (w(ue ) )2 ≤ n−1 . ue )ij = ωe Se ∈Fe

Similarly, for each half-edge h ∈ H, let Fh be either S1 or Sm+1 , and let ch be the corresponding coefficient vector, either b1 or b2 . The endpoint factor attached to h has the form X gh (i) = ηh ch (Sh )xSi h , ηh ≤ n−1/2 . Sh ∈Fh

For each admissible assignment of distinct sample labels to the r blocks, independence across the distinct block labels factors the expectation over the vertices of G♯ . Since the number of such assignments is at most nr , taking absolute values gives the reduction Y Y E[Sπ♯ ] ≤ |Θ| nr ωe ηh Γ, (212) e∈E

h∈H

where the remaining feature sum Γ is defined by Y X Y Γ := |ch (Sh )| µT (Sf )f ∈∂T . T ∈π ♯

(Sf )f ∈E∪H h∈H

Here, for each block T ∈ π ♯ , the local Walsh moment is defined as hY i µT (Sf )f ∈∂T := E xSf . f ∈∂T

Since each µT is a Fourier–Walsh moment, it satisfies the pointwise bound 0 ≤ µT (Sf )f ∈∂T ≤ 1. We next bound Γ. Since all feature families are subcollections of 2[d] , and since the number of surviving edges and half-edges depends only on m, the internal feature sums contribute only an Od (1) factor. This gives the estimate Y X Γ ≤ Od (1) |ch (Sh )|. h∈H Sh ∈Fh

For each endpoint feature family, Cauchy–Schwarz gives the bound X |ch (Sh )| ≤ |Fh |1/2 ∥ch ∥2 ≤ Od (∥ch ∥2 ). Sh ∈Fh

Combining the preceding two estimates yields ! Γ ≤ Od

∥ch ∥2

(213)

h∈H

It remains to account for the endpoint coefficients. In the original second-moment expansion, there are exactly two b1 -endpoint factors and two b2 -endpoint factors. Each such factor either remains as a half-edge or has already been absorbed into the external scalar Θ. This accounting gives the endpoint bound Y |Θ| ∥ch ∥2 ≤ ∥b1 ∥22 ∥b2 ∥22 . h∈H

Together with (213), this bound implies |Θ| Γ ≤ Od ∥b1 ∥22 ∥b2 ∥22 .

(214)

We now count the powers of n. Each surviving sample position appears in exactly two surviving factors. The corresponding degree count is the identity 2|E| + |H| = 2q. Since ωe ≤ n−1 for every internal edge and ηh ≤ n−1/2 for every half-edge, the product of weights satisfies Y Y ωe ηh ≤ n−|E|−|H|/2 = n−q . (215) e∈E

h∈H

Substituting (214) and (215) into (212) gives E[Sπ♯ ] ≤ nr−q Od ∥b1 ∥22 ∥b2 ∥22 . The counting inequality (211), together with r ≥ 1, implies r − q ≤ −r ≤ −1. The preceding estimate therefore gives E[Sπ♯ ] = Od (n−1 )∥b1 ∥22 ∥b2 ∥22 . This proves the proposition. 92

Technical lemmas and propositions

Proposition 43. Let X be a random element of Hn×d , and let f1 , . . . , fd : Hn×d → R. Assume hypercontractivity: for every p ≥ 2 and every i ∈ [d], ∥fi (X)∥Lp ≤ Cp,ℓ ∥fi (X)∥L2 . Then, for every ϵ > 0, there exists a constant Cϵ,ℓ , depending only on ϵ and ℓ, such that 2

max fi (X) i∈[d]

≤ Cϵ,ℓ dϵ max ∥fi (X)∥2L2 . i∈[d]

2 More explicitly, one may take Cϵ,ℓ = C2q,ℓ for any choice of q > 1 satisfying 1/q ≤ ϵ.

Proof of Proposition 43. Fix ϵ > 0, and choose q > 1 such that 1/q ≤ ϵ; for instance, one may take q = max{2, 1/ϵ}. Define M (X) := max fi (X). i∈[d]

For any real numbers a1 , . . . , ad , the inequality (maxi ai )2 ≤ maxi a2i holds. Applying this pointwise gives the initial reduction h i ∥M ∥2L2 = E M (X)2 ≤ E max fi (X)2 . (216) i∈[d]

Set

Y (X) := max fi (X)2 . i∈[d]

The random variable Y is nonnegative. By Lyapunov’s inequality, or equivalently the monotonicity of Lp -norms on a probability space, the L1 -norm of Y is bounded by its Lq -norm. Thus the right-hand side of (216) satisfies h i 1/q E max fi (X)2 = ∥Y ∥L1 ≤ ∥Y ∥Lq = E[Y (X)q ] . i∈[d]

Using the definition of Y , this bound can be written as E[Y (X)q ] The elementary inequality maxi bi ≤ bound

1/q

h i1/q = E max |fi (X)|2q . i∈[d]

i=1 bi , valid for nonnegative numbers bi , gives the moment

d h i1/q h X i1/q E max |fi (X)|2q ≤ E |fi (X)|2q i∈[d]

i=1 d X 1/q = E |fi (X)|2q i=1

1/q ≤ d1/q max E |fi (X)|2q . i∈[d]

(217)

For each i ∈ [d], the assumed hypercontractive estimate, applied with exponent 2q, gives ∥fi (X)∥L2q ≤ C2q,ℓ ∥fi (X)∥L2 . Equivalently, after squaring both sides, this estimate becomes 1/q 2 E |fi (X)|2q = ∥fi (X)∥2L2q ≤ C2q,ℓ ∥fi (X)∥2L2 .

(218)

Combining (216), (217), and (218) gives the estimate 2

max fi (X) i∈[d]

2 ≤ d1/q C2q,ℓ max ∥fi (X)∥2L2 . i∈[d]

2 Since 1/q ≤ ϵ and d ≥ 1, the factor d1/q is bounded above by dϵ . Therefore, with Cϵ,ℓ := C2q,ℓ , the desired estimate follows.

Proposition 44 (Weighted sums of Fourier–Walsh polynomials). Let 1 ≤ m ≤ q, let p = (p1 , . . . , pq ) ∈ Nq , and fix S ∗ ⊂ [d] with |S ∗ | = p∗ . For each i ∈ [q], let Si be a collection of subsets of [d] such that |S| ≤ pi for every S ∈ Si . Then, for every choice of vectors a[i] ∈ R|Si | , i = 1, . . . , m, one has X

m h i Pq Y ∗ [1] [m] aS1 · · · aSm E xS1 · · · xSq xS ≤ Cq,p d j=m+1 pj /2 ∥a[i] ∥2 .

S1 ∈S1 ,...,Sq ∈Sq

(219)

i=1

Proof of Proposition 44. Let µ denote the underlying probability measure, and write ∥ · ∥Lr for the corresponding Lr (µ)-norm. For each i ∈ [m], define the weighted Walsh polynomial fi by X [i] fi (x) := aS xS . S∈Si

For each j = m + 1, . . . , q, define the unweighted Walsh polynomial gj by X gj (x) := xS . S∈Sj

Expanding these definitions shows that the expression inside the absolute value in (219) is the quantity   q m Y Y ∗ L := Eµ xS fi (x) gj (x) . (220) i=1

j=m+1

∗

The identity |xS | = 1 holds pointwise, so multiplication by xS preserves every Lr (µ)-norm. An application of Cauchy–Schwarz gives the estimate m Y

|L| ≤

i=1

To estimate this L2 -norm, set

( hr :=

q Y

j=m+1

. L2

fr , 1 ≤ r ≤ m, gr , m + 1 ≤ r ≤ q.

Hölder’s inequality, applied with exponent q to the q factors |hr |2 , yields the bound q Y r=1

Eµ

q Y

!1/2 2

|hr |

≤

r=1

q Y

∥hr ∥L2q .

r=1

Combining the preceding two estimates gives the key reduction |L| ≤

m Y i=1

q Y

∥fi ∥L2q

j=m+1

∥gj ∥L2q .

(221)

It remains to bound the factors on the right-hand side of (221). Since fi is a Walsh polynomial of degree at most pi , the Bonami hypercontractive inequality gives the estimate ∥fi ∥L2q ≤ (2q − 1)pi /2 ∥fi ∥L2 . The orthonormality of the Fourier–Walsh basis gives the identity X [i] ∥fi ∥2L2 = |aS |2 = ∥a[i] ∥22 . S∈Si

Combining these two estimates yields ∥fi ∥L2q ≤ (2q − 1)pi /2 ∥a[i] ∥2 .

(222)

The same argument applies to the unweighted factors. Since gj is a Walsh polynomial of degree at most pj , the Bonami hypercontractive inequality gives the estimate ∥gj ∥L2q ≤ (2q − 1)pj /2 ∥gj ∥L2 . The orthonormality of the Fourier–Walsh basis gives the identity X ∥gj ∥2L2 = 1 = |Sj |. S∈Sj

Since every set in Sj has cardinality at most pj , the cardinality of Sj satisfies the bound |Sj | ≤

pj X d r=0

≤ Cpj dpj .

Combining the preceding three estimates gives ∥gj ∥L2q ≤ Cq,p dpj /2 . Substituting (222) and (223) into (221) yields the final estimate ! m Pq Y [i] |L| ≤ Cq,p ∥a ∥2 d j=m+1 pj /2 . i=1

By the representation of L in (220), this estimate is exactly the desired bound (219).

(223)

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