¹

1 Introduction↩︎

The debiased machine learning (DML), also known as the double/debiased machine learning framework, has become a leading approach to two-step estimation with high-dimensional or nonparametric nuisance components; see, among many others, [1], [2], and [3]. A central feature of this literature is cross-fitting, whereby sample splitting is used to separate nuisance estimation from target parameter evaluation. This device is motivated by two considerations. First, the conventional view holds that it mitigates overfitting bias arising from the use of highly flexible first-stage learners. Second, it relaxes empirical process conditions by reducing the dependence between first-stage estimation errors and second-stage score evaluation. As a result, cross-fitting has become close to a default recommendation in both theoretical analyses and empirical implementations of DML procedures.

At the same time, many empirical applications in economics and finance involve multiway clustered dependence; see, for example, [4], [5]. In such settings, observations may be correlated along multiple dimensions, such as firm and time or region and industry. Extensions of DML to multiway clustered environments are developed in [6].

However, extensive sample splitting is not without cost. First, with a finite number of folds, cross-fitting yields a random estimator, which may hinder reproducibility. Increasing the number of folds can mitigate this concern to some extent, but at the expense of substantially greater computational burden, particularly when the first stage is complex and/or requires extensive tuning. Second, cross-fitting effectively reduces the sample size available to each first-stage problem. This is especially consequential when these stages involve high-dimensional or nonparametric estimation, where performance is inherently variance-sensitive. The resulting loss in nuisance estimation precision may, in finite samples, translate into non-negligible efficiency losses for the parameter of interest.

These concerns are further amplified under multiway clustering, where the effective sample size is determined by the number of independent cluster units rather than the total number of observations. Partitioning the data into folds may therefore leave each subsample with only a limited number of independent clusters, making cross-fitting particularly costly, in addition to increasing computational burden. Moreover, in two-step procedures such as double machine learning, overfitting bias arising from first-stage nuisance estimation need not be intrinsically detrimental for inference on the target parameter, in contrast to classical one-step settings. This suggests that the conventional rationale for cross-fitting may be less central than is sometimes presumed. Indeed, even under i.i.d. settings, the literature provides little general theoretical support for the advantages of cross-fitting, apart from a few special cases considered, for example, by [7].

Motivated by these considerations, a growing literature seeks to weaken or dispense with cross-fitting in general nonlinear estimation problems. One approach to obtaining theoretical guarantees for DML without cross-fitting is the “localisation" method, as employed for example in [8], [9]. The key idea is to analyse the supremum of an empirical process indexed by estimating equations evaluated over deterministic sets that localise the nuisance parameter around its true value. These sets are constructed so that the nuisance estimator lies in them with probability approaching one, thereby replacing a stochastic index with a deterministic one and disentangling the dependence between estimating equations and first-stage estimates. If the localisation sets shrink at appropriate rates—typically verified through maximal inequalities—the associated empirical process remainder is of smaller order than the leading asymptotically linear term and does not affect the limiting distribution. An alternative strategy is based on a”stability" conditions; see [10] under i.i.d. and a generalisation in [11] under spatial/network \(\beta\)-mixing. Verifying such conditions, however, is often substantially more involved beyond certain well-understood cases, and we therefore do not pursue this route.

In this paper, we develop a general asymptotic theory for two-step debiased GMM estimators under multiway clustered dependence, without relying on cross-fitting. We consider a broad class of locally robust two-step GMM problems similar to those studied in [2] in which low-dimensional target parameters are identified by orthogonal moment conditions that depend on high-dimensional or nonparametric nuisance components. The data are allowed to exhibit dependence along an arbitrary number of clustering dimensions, accommodating empirical settings in which correlation arises simultaneously across, for example, firms, time periods, locations, or networks. Our results establish asymptotic linearity and asymptotic normality of the resulting estimators under conditions that permit flexible highly first-stage learners while avoiding sample splitting. The analysis explicitly accounts for the reduced effective sample size induced by multiway clustering and provides inference procedures that remain valid when the number of independent cluster units, rather than the total number of observations, governs the stochastic order. By combining orthogonality with a localisation-based empirical process argument tailored to multiway clustered arrays, we show that the impact of first-stage estimation can be controlled without cross-fitting. This yields a unified framework for debiased GMM inference that is well suited to empirically relevant clustered environments where conventional cross-fitting can be statistically and computationally costly.

Maximal inequalities are central to localisation-based arguments in DML, yet for multiway clustered—more precisely, separately exchangeable (SE)—arrays, the available theory remains limited. In particular, no general global maximal inequality accommodates arbitrary numbers of clustering dimensions \(K\), arbitrary moments \(q\), and infinite pointwise measurable function classes. Existing results address only special cases: Theorem B.2 of [12] allows general \(K\) and \(q\in[1,\infty)\) but restricts attention to finite classes, while Lemma C.3 of [13] covers general classes only for \(q=1\) and \(K=2\). The situation is even more restrictive for local maximal inequalities, which are essential for sharper convergence rates: unlike in the related \(U\)-statistics literature (see [14]), beyond the \(K=1\) (i.i.d.) case, no such result appears to be available for SE arrays. These gaps reflect the intrinsic difficulty posed by multiway dependence, which generates complex interactions across observations and undermines classical tools such as Hoeffding averaging. To overcome this, we develop a new proof strategy based on a transversal partition of the index set that effectively decouples dependence and permits the use of the Hoffmann–Jørgensen inequality, yielding sharp higher-moment bounds. This approach delivers both global and local maximal inequalities for potentially uncountable, pointwise measurable function classes under SE sampling.

The paper is organised as follows. Section 2 introduces the multiway clustered sampling setups and notations. Section 3 develops the cross-fitting-free debiased GMM estimator and presents the main asymptotic results under general high-level conditions, while also providing three examples for the rate and complexity conditions and validity of variance estimation. Section 4 provides the core technical tools in the form of new maximal inequalities for empirical processes for separately exchangeable arrays. Proofs and supplementary arguments are collected in the appendices.

2 Setups and Notations↩︎

In this section, we introduce the framework of multiway clustered sums that will be used throughout the paper. Let \(K\) be a fixed positive integer, and denote a \(K\)-tuple index by \(\boldsymbol{i} = (i_1, i_2, \dots, i_K) \in \mathbb{N}^K.\) Suppose we observe a \(K\)-way array of data \[\{ X_i : \boldsymbol{i} \in [N_1] \times \cdots \times [N_K] \},\] where \([N_k] := \{1,\dots,N_k\}\) denotes the index set of dimension \(k\) and \(N_k\) is the sample size. Let \(\boldsymbol{N} = (N_1, N_2, \dots, N_K)\) and define \([\boldsymbol{N}] = \prod_{k=1}^{K} \{1,2,\dots,N_k\}.\) We also denote \[N = \prod_{k=1}^K N_k,\quad n = \min\{N_1, N_2, \dots, N_K\},\quad \text{and}\quad \overline{N} = \max\{N_1,...,N_K\}.\]

Suppose \(\{ X_i\}_i\) are random variables defined on a probability space \((S,\mathcal{S}, P)\)², and \(\{ X_i\}\) satisfy the separate exchangeability (SE) and dissociation (D) conditions defined below.

1. For any \(\pi =(\pi_1,\dots,\pi_K)\), a \(K\)-tuple of permutations of \(N\), \(\{X_{\boldsymbol{i}}\}_{\boldsymbol{i} \in [\boldsymbol{N}]}\) and \(\{X_{\pi(\boldsymbol{i})}\}_{\boldsymbol{i} \in [\boldsymbol{N}]}\) are identically distributed.

2. For any two disjoint sets of indices \(I,I'\subset \mathbb{N}^K\), \(\{X_{\boldsymbol{i}}\}_{{\boldsymbol{i}}\in I}\) and \(\{X_{\boldsymbol{i}}\}_{{\boldsymbol{i}}\in I'}\) are independent.

Under Conditions (SE) and (D), the Aldous–Hoover–Kallenberg (AHK) representation (see Corollary 7.35 in [15]) guarantees the existence of the following representation: \[\begin{align} X_{\boldsymbol{i}} = \tau\Bigl( \{ U_{\boldsymbol{i} \odot \boldsymbol{e}} \}_{\boldsymbol{e} \in \{0,1\}^K \setminus \{\boldsymbol{0}\}} \Bigr), \label{eq:AHK95representation} \end{align}\tag{1}\] where \(\odot\) denotes the Hadamard (element-wise) product, the collection \(\{ U_{\boldsymbol{i} \odot \boldsymbol{e}} : \boldsymbol{i} \in \mathbb{N}^K,\; \boldsymbol{e} \in \{0,1\}^K \setminus \{\boldsymbol{0}\} \}\) consists of mutually independent and identically distributed (i.i.d.) random variables, and \(\tau\) is a Borel measurable map taking values in \(\mathcal{S}\).

We say a class of functions \(\mathcal{F}:\mathcal{S}\to \mathbb{R}\) is pointwise measurable if there exists a countable subclass \(\mathcal{F}'\subset \mathcal{F}\) such that for each \(f\in \mathcal{F}\), there exists a sequence \((f_j)_j\subset \mathcal{F}'\) such that \(f_j\to f\) pointwisely. Given the observed set of random variables \(\{X_{\boldsymbol{i}}:{\boldsymbol{i}}\in [\boldsymbol{N}]\}\) that satisfy Conditions (SE) and (D), and a pointwise measurable class of functions \(\mathcal{F}\) with elements \(f: \mathcal{S}\to \mathbb{R}\), define the sample mean process by \(\mathbb{E}_N f = N^{-1}\sum_{\boldsymbol{i} \in [\boldsymbol{N}]} f(X_{\boldsymbol{i}})\) and, suppose \(P|f|:=\int |f| dP<\infty\)³, the empirical process by \[\mathbb{G}_{n}(f) = \frac{\sqrt{n}}{N} \sum_{\boldsymbol{i} \in [\boldsymbol{N}]} \Bigl\{ f(X_{\boldsymbol{i}}) - P(f)\bigr] \Bigr\}.\] Throughout, we use the shorthand \(\mathbb{E}[f(X_{\boldsymbol{1}})]= P(f)\), where \(\boldsymbol{1}=(1,...,1)\), for expectations taken with respect to the data-generating distribution.

Notation.↩︎

Let \(\mathbb{N}\) denote the set of positive integers and \(\mathbb{R}\) for the real line. For \(a,b\in \mathbb{R}\), let \(a\vee b=\max\{a,b\}\) and \(a\wedge b = \min\{a,b\}\). Denote for \(m\in \mathbb{N}\) that \([m] = \{1,2,\ldots,m\}.\) For real vectors \(\boldsymbol{a}= (a_{1},\dots,a_{K})\) and \(\boldsymbol{b} = (b_{1},\dots,b_{K})\), we denote \(\boldsymbol{a} \le \boldsymbol{b}\) for \(a_{j} \le b_{j}\) for all \(1 \le j \le K\). Let \(\mathrm{supp}(\boldsymbol{a}) = \{ j : a_j \ne 0\}\). We denote by \(\odot\) the Hadamard product: for \({\boldsymbol{i}}= (i_1,\dots,i_K)\) and \(\boldsymbol{j} = (j_1,\dots,j_K)\), \({\boldsymbol{i}}\odot \boldsymbol{j}= (i_1 j_1,\dots,i_K j_K)\). For each \(k=0,1,2,...,K\), define \(\mathcal{E}_k=\{{\boldsymbol{e}}\in \{0,1\}^K: \Vert{\boldsymbol{e}}\Vert_0 =k\}\) and thus \(\{0,1\}^K=\cup_{k=0}^K \mathcal{E}_k\). For \(q\in [1,\infty]\), let \(\|f\|_{Q,q}=( Q|f|^q)^{1/q}\). For a non-empty set \(T\) and \(f:T\to \mathbb{R}\), denote \(\|f\|_T=\sup_{t\in T}|f(t)|\). For a pseudometric space \((T,d)\), let \(N(T,d,\varepsilon)\) denote the \(\varepsilon\)-covering number for \((T,d)\). We say \(F:\mathcal{S}\to \mathbb{R}_+\) is an envelope for a class of functions \(\mathcal{F}\ni f:\mathcal{S}\to \mathbb{R}\) if \(\sup_{f\in\mathcal{F}}|f(x)|\le F(x)\) for all \(x\in\mathcal{S}\).

3 Debiased Machine Learning for GMM↩︎

In this section, we DML type two-step estimation and inference approaches in a setting where data is multiway clustered. Particularly, it is of practical importance to study debiased machine learning without sample-splitting due to the poor usage of samples when splitting the data. While cross-fitting can improve the sample usage by switching the roles of split samples, the improvement is limited in a multiway clustered setting where cross-fitting is done in a way that more data is excluded when estimating the high-dimensional nuisance parameters.

Since most econometric and statistical models can be reduced to moment restrictions, we consider DML without sample splitting in a GMM setup similar to those considered in [2]:

1. The parameters \(\theta_0 \in \Theta \subset R^d\) of interest are fixed-dimensional and (over-)identified by a set of moment conditions that depends on high-dimensional nuisance parameters: \(g: \mathcal{X}\times\Theta\times {\Gamma}\rightarrow\mathbb{R}^q\) such that \[\label{gmm95moment} \mathbb{E}\left[ g(X,\theta_0,\gamma_0) \right] = 0,\tag{2}\]

2. The nuisance parameters \(\gamma_0\in \Gamma\) are exactly identified and estimated by machine learners using the full sample.

3. With the full sample again, solving the GMM problem using an orthogonalized moment function and the corresponding optimal weighting matrix, with plug-in nuisance estimates.

DML for two-step GMM with i.i.d. data is studied in [2]. In contrast, our setting involves multiway-clustered sampling, which leads to substantially different asymptotic arguments. Moreover, our theoretical framework eliminates the need for cross-fitting.

3.1 Main results for the debiased GMM estimator.↩︎

The first component of DML is the orthogonalisation of the moment condition. Specifically, we construct \(\psi(X,\theta_0,\eta_0)\) by adding an adjustment term (which may depends on extra nuisance parameters contained in \(\eta_0\in \Gamma\)) to \(g\) such that \(\psi\) is mean zero and the path-wise derivative with respect to \(\eta\) in the direction \(\widetilde{\eta}\in \Gamma\) is zero (or vanishing) when evaluated at the truth: \[\begin{align} \partial_\eta \mathbb{E}[\psi(X,\theta_0,\eta_0)](\widetilde{\eta}) := \partial_\tau \mathbb{E}[\psi(X,\theta_0,\eta_0 + \tau\widetilde{\eta})]_{\tau = 0} = 0. \label{orthogonal} \end{align}\tag{3}\] Such adjustment offsets the effect of local perturbation of \(\gamma\) on the identifying moment condition. This component of DML is a property with respect to the population moment condition, i.e., irrelevant of multiway clustering or cross-fitting, and it is well-established in the GMM setting due to aforementioned literature, among others. Therefore, we take this condition as given for our analyses.

Let \(\widehat\eta\) be some machine learners that are appropriate for multiway clustering data⁴, and let \(\widehat\psi_N\) denote the empirical average with plug-in estimate \(\widehat\eta\): \[\begin{align} \widehat\psi_N(\theta) = \mathbb{E}_N[\psi(X,\theta,\widehat\eta)] \end{align}\] With some positive semi-definite weighting matrix \(\widehat\Upsilon\) (e.g., the inverse of a multiway cluster-robust variance-covariance estimator of \(\psi(X;\widehat\theta^{(0)},\widehat\eta)\) with some initial estimate \(\widehat\theta^{(0)}\)), the debiased GMM estimator of \(\theta_0\) is defined as \[\begin{align} \widehat\theta = \arg \min_{\theta\in \Theta} \widehat\psi_N'(\theta) \widehat\Upsilon \widehat\psi_N(\theta). \label{gmm95estimator} \end{align}\tag{4}\] When \(\theta_0\) is exactly identified, the debiased GMM estimator reduces to \(\widehat\theta\) as a solution to \[\begin{align} \widehat\psi_N(\theta) = 0. \end{align}\] Our analyses are based on the general case (4 ).

We denote the population and empirical Jacobian as \[J_0 := -\partial_\theta \mathbb{E}[\psi(X;\theta,\eta_0)]\big|_{\theta=\theta_0}, \qquad \widehat J_N(\widehat\theta) := -\partial_\theta \widehat\psi_N(\theta)\big|_{\theta=\widehat\theta}.\] Suppose we have an interior minimizer \(\widehat\theta\) from (4 ), then the first-order condition holds as follows: \[\begin{align} 0=\widehat J_N(\widehat\theta)' \widehat\Upsilon \widehat\psi_N(\widehat\theta) \label{foc} \end{align}\tag{5}\]

Let \(f(\eta) = \psi(X,\theta_0,{\eta}) - \psi(X,\theta_0,{\eta}_0)\). By a standard mean-value expansion of \(\psi_N(\widehat\theta)\) in (5 ), we can write \(\sqrt{n}(\widehat\theta - \theta_0)\) as a function of a well-behaved term \(\mathbb{E}_N [\psi(X,\theta_0,{\eta}_0)]\), an empirical process term \(\mathbb{G}_{n}\left(f(\widehat\eta)\right)\), an extra error term \(\mathbb{E}[ f(\eta)]_{\eta = \widehat\eta}\) due to the nuisance parameter estimation, as well as the empirical Jacobian \(J_N(\widehat\theta)\) and feasible weighting matrix \(\widehat\Upsilon\). As in the DML literature, \(\mathbb{E}[ f(\eta)]_{\eta = \widehat\eta}\) can be bounded by the orthogonality condition. It is relatively straightforward to deal with the empirical Jacobian term, and it is standard in the literature to put aside the weighting matrix estimation, as long as \(\widehat{\Upsilon}\overset{p}{\to}\Upsilon\) to some positive-definite matrix \(\Upsilon\). Now the difficult term left is \(\mathbb{G}_{n}(f(\widehat\eta))\), and we bound it using a localisation approach through the maximum inequality under multiway clustering.

The idea of the localisation approach is that under the exchangeability and dissociation conditions, we can utilize the Hoeffding decomposition of the empirical process \(\mathbb{G}_{n}\left(f(\widehat\eta)\right)\) and bound each of the decomposed terms by the maximum inequality in a neighborhood of \(\eta_0\). As long as the first-step machine learner of the nuisance parameter lies in the neighborhood \(\Gamma_n(\eta_0)\) with high probability, we can show \(\mathbb{G}_{n}\left(f(\widehat\eta)\right)\) vanishes asymptotically. To formally define the neighborhood of \(\eta_0\in\Gamma\), we equip \(\Gamma\) with the \(L_2(P)\) norm \(\Vert.\Vert_{P,2}\).

Assumption 1(i) characterizes the multiway clustered data by the exchangeability and dissociation conditions, which are standard in clustering robust inference literature. Assumptions 1(ii) and (iii) are score regularity conditions. Assumption 1(ii) is satisfied when the scores come from a likelihood function or moment conditions that are smooth in terms of the nuisance parameters. Assumption 1(iii) is a nonlinear counterpart of the linear-in-\(\theta\) condition common in the DML literature. For a score that is twice differentiable in \(\theta\), the existence of the integrable envelope \(B(X,\eta)\) reduces to the integrability of the Hessian matrix locally. See Remark 1 below for more details on the choice of \(B(X,\eta)\). The locality in \(\Gamma_n(\eta_0)\) ensures that \(\eta\) takes values that do not explode up the envelope, e.g., \(\eta\) as inverse probability weights. Assumption 1(iv) is a high-level condition governing the quality of nuisance parameter estimation. This requirement can be verified using existing theoretical results for a range of machine-learning estimators under multiway clustering; for instance, in the case of LASSO, it follows from Proposition 2 of [12]. Assumptions 1(v) and (vi) are standard and mild finite moment and full rank conditions.

Following Chapter 3.6 in [19], a function class \(\mathcal{F}\) on \(\mathcal{S}\) with a measurable envelope \(F\) is called Vapnik–Chervonenkis-type (VC-type) with characteristics \((A,v)\) if \[\begin{align} \sup_Q N(\mathcal{F},\|\cdot\|_{Q,2},\varepsilon\|F\|_{Q,2})\le \left(\frac{A}{\varepsilon}\right)^v\: \text{ for all } 0<\varepsilon\le 1,\label{eq:VC-type} \end{align}\tag{6}\] where the supremum is taken over all finite discrete distributions. It can be shown that a wide range of commonly used models and estimators in econometrics, machine learning, and statistics give rise to sequences of function classes that satisfy this VC-type condition; see Section 3.2 below for illustrative examples.

The following theorem presents our first main result, establishing the asymptotic linearity and asymptotic normality of a generic debiased GMM estimator without cross-fitting.

A proof can be found in Section 6.1 in the appendix. The additional condition (b) in statement (2) of the theorem is a non-degeneracy requirement, ensuring that at least one clustering dimension enters the score in a linear manner. This condition is mild for larger \(K\)’s, as it only requires that at least a single one latent shock of the \(K\) clustering dimensions has a non-trivial effect on \(\psi\). This condition was also imposed in, e.g., [20], [6] and [12]. In the case of i.i.d data, the non-degeneracy condition does not hold. In such conventional settings, the asymptotic normality result for the full-sample DML approach has been established in [8], among others, while here we focus on the non-degenerate case.

To build intuition, note that the first-order condition implies the expansion \[\begin{align} \sqrt{n}M^{-1}(\widehat\theta - \theta_0) \;\approx\; \mathbb{E}_N[\psi(X,\theta_0,\eta_0)] \,+\, \sqrt{n}\,\mathbb{E}[f(\eta)]_{\eta=\widehat\eta}\, +\, \mathbb{G}_n\bigl(f(\widehat\eta)\bigr), \end{align}\] for some invertible matrix \(M\), where \(f(\eta)=\psi(X,\theta_0,\eta)-\psi(X,\theta_0,\eta_0).\) Such a decomposition is standard in semiparametric theory; see, for example, [21].

The first term is asymptotically normal. The second term is controlled by the orthogonality condition 3 , and is therefore negligible. Consequently, the main technical challenge is to control the localised empirical process \(\mathbb{G}_n\bigl(f(\widehat\eta)\bigr),\) which is non-standard due to the dependence of \(\widehat\eta\) on the full sample. A standard approach is cross-fitting, which removes this dependence. Conditional on \(\widehat\eta\), one may apply Hoeffding-type decomposition and Markov’s inequality to obtain, under suitable smoothness conditions, with probability \(1-o(1)\) \[\mathbb{G}_n\bigl(f(\widehat\eta)\bigr) \;\lesssim\; \|\widehat\eta - \eta_0\|^{v}_{P,2} \quad \text{for some } v>0.\]

An alternative is a localisation argument. Suppose there exists a sequence of shrinking function classes \(\{\mathcal{F}_n\}\) such that \(\mathbb{P}(\widehat\eta \in \mathcal{F}_n) \to 1.\) Then, with probability \(1-o(1)\), \[\bigl|\mathbb{G}_n(f(\widehat\eta))\bigr| \;\le\; \sup_{\eta\in\mathcal{F}_n} \bigl|\mathbb{G}_n(f(\eta))\bigr|,\] which removes the stochastic dependence on \(\widehat\eta\). Such classes \(\{\mathcal{F}_n\}\) can often be constructed tightly when the convergence rate of \(\|\widehat\eta - \eta_0\|_{P,2}\) is available. This is typically the case, as the same rate is also needed to verify the orthogonality condition regardless of whether cross-fitting is used.

In the classical semiparametric literature, the function class \(\mathcal{F}\) is typically fixed, and stochastic equicontinuity follows from standard uniform (functional) CLT arguments. In contrast, with machine-learning first stages, a fixed \(\mathcal{F}\) is generally too large to control, necessitating shrinking (localised) classes. This localisation strategy underlies the i.i.d.analyses of [8] further generalised in [9], which rely on maximal inequalities from [22] to control the supremum. Since comparable results are unavailable under multiway clustering, we develop the required global and local maximal inequalities in Section 4.

3.2 Verification of complexity and rate conditions: three examples↩︎

Under Assumption 1, to apply Theorem 1 it suffices to verify the high-level VC-type condition in 6 and the rate condition in 8 . Verifying these conditions is not entirely straightforward in general, owing to their abstract nature. Below we discuss several examples of machine learning estimators for the first-stage nuisance function \(\eta\) that satisfy these two conditions.

To isolate the role of the first-stage nuisance parameter learner, we consider \(\psi(\cdot,\eta)\) as a map in \(\eta\) that preserves the VC-type properties of the underlying function class \(\mathcal{G}_n\) to which \(\eta\) belongs. For instance, \(\psi(\cdot,\eta)\) may be a monotone or Lipschitz transformation of \(\eta\), or a finite combination such as sums, products, minima, or maxima; see Section 3.6 of [19].

A proof can be found in Section 6.2 in the appendix.

Case (i) admits generalized linear sparse models, including sparse linear, logit, and exponential models as special cases. These models correspond to LASSO-type (\(\ell^1\)-penalty) machine learners for a sparse generalised linear model. Cases (ii) and (iii) correspond to the regression tree and deep neural networks, respectively. More details are given in the appendix on how the rate conditions are obtained. Basic definitions and textbook treatments of these methods can be found in e.g. [26].

3.3 Variance estimation.↩︎

For hypothesis testing using the results given in Theorem 1, a missing piece is the unknown asymptotic variance \(V\). In this section, we propose a full-sample variance estimator that takes into account (1) multiway clustering dependence and (2) estimation errors from both high-dimensional nuisance estimation and the GMM estimation. To account for the multiway clustering dependence, we follow the formulation of the multiway clustering-robust variance estimator in [27], except that the empirical scores here are replaced by the Neyman orthogonalised scores with estimated nuisance parameters.

For any \(\boldsymbol{i}, \boldsymbol{j} \in [\boldsymbol{N}]\), let \(i_k\) and \(j_k\) denote their \(k\)-th elements, and let \(\mathbb{1}_{k}\{\boldsymbol{i}, \boldsymbol{j}\}\) indicate whether the two observations \(\boldsymbol{i}, \boldsymbol{j}\) share the same cluster at \(k\)-th dimension, i.e., \(\mathbb{1}_{k}\{\boldsymbol{i}, \boldsymbol{j}\} = \mathbb{1}\{ i_k = j_k \}.\) We define the estimator for the middle term \(\Psi_0\) as follows: \[\begin{align} \widehat{\Psi}_N(\widehat\theta) = & \sum_{ k=1 }^K \widehat\Psi_{N,k}(\widehat\theta), \\ \widehat\Psi_{N,k}(\widehat\theta)= & \frac{n}{N^2} \sum_{\boldsymbol{i,j\in [N]}} \psi(X_{\boldsymbol{i}},\widehat\theta,\widehat\eta) \psi(X_{\boldsymbol{j}},\widehat\theta,\widehat\eta)' \mathbb{1}_{k}\{\boldsymbol{i}, \boldsymbol{j}\} \end{align}\] Then the estimator for \(V\) is given as follows: \[\begin{align} \label{eq:var95est} \widehat{V}=\left(\widehat{J}_N(\widehat\theta)' \widehat{\Upsilon} \widehat{J}_N(\widehat\theta)\right)^{-1} \widehat{J}_N(\widehat\theta)' \widehat{\Upsilon} \widehat{\Psi}_N(\widehat\theta) \widehat{\Upsilon} \widehat{J}_N(\widehat\theta) \left(\widehat{J}_N(\widehat\theta)' \widehat{\Upsilon} \widehat{J}_N(\widehat\theta)\right)^{-1}. \end{align}\tag{9}\]

In practice, if there are more than one observation in some cells \(\boldsymbol{i}\), we simply aggregate within each cell by replacing \(\psi(X_{\boldsymbol{i}},\widehat{\theta},\widehat{\eta})\) with the sum of empirical scores \(\psi\) within that cell. Since this generalisation would not change the main analysis except for complications in notations, we focus on the case with exactly one observation in each cell.

A proof can be found in Section 6.3 in the appendix.

Theorem 3 establishes the consistency of the variance estimator using the full sample, under the same non-degeneracy condition as in the second statement of Theorem 1. The extra moment conditions mildly strengthen the moment conditions in Theorem 1. As in Theorem 1, these local integrability conditions can be delivered by integrability conditions of the score, Jacobian, and the Hessian, given enough smoothness in \(\eta\).

\(\widehat{V}\) is positive semi-definite by construction because each \(\widehat\Psi_{N,k}(\widehat\theta)\) is positive semi-definite mechanically. This can be seen easily in the case \(K=2\), in which case each \(\widehat\Psi_{N,k}(\widehat\theta)\) reduces to a one-way cluster variance estimator. A caveat is that when none of the cluster matters, e.g., i.i.d. data, the non-degeneracy condition can fail, and this variance estimator would overestimate the asymptotic variance \(V\) and result in a conservative test, which is well-known in the cluster robust inference literature (e.g., see [28]). A potential fix for this issue is to remove double-counting terms in \(\widehat\Psi_{N}(\widehat\theta)\) by defining \[\begin{align} \widetilde{\Psi}_N(\widehat\theta) = & \sum_{ \boldsymbol{e}\in\{0,1\}^K, \Vert \boldsymbol{e}\Vert=r } (-1)^{r+1}\widetilde{\Psi}_{N,\boldsymbol{e}}(\widehat\theta), \\ \widetilde{\Psi}_{N,\boldsymbol{e}}(\widehat\theta) =& \frac{n}{N^2} \sum_{\boldsymbol{i,j\in [N]}} \psi(X_{\boldsymbol{i}},\widehat\theta,\widehat\eta) \psi(X_{\boldsymbol{j}},\widehat\theta,\widehat\eta)' I_{\boldsymbol{e}}\{\boldsymbol{i}, \boldsymbol{j}\}, \end{align}\] where \(I_{\boldsymbol{e}}\{\boldsymbol{i}, \boldsymbol{j}\}= \mathbb{1}\{\boldsymbol{e}\odot \boldsymbol{i} = \boldsymbol{e}\odot \boldsymbol{j}\}\), which instead indicates whether the two observations share the same clusters over the whole support of \(\boldsymbol{e}\). This is basically a DML version of the \(K\)-way generalisation of variance estimator proposed in [29], referred to as the CGM estimator⁷. In a parametric setting, [27] shows that these two types of variance estimators are both consistent for the asymptotic variance under non-degeneracy. However, the CGM estimator involves more terms to calculate and is not guaranteed to be positive semi-definite. In practice, it is rarely true that multi-dimensional data is i.i.d because of the common existence of unobserved heterogeneous effects.

For some degenerate yet cluster-dependent scenarios, such as those studied by [31], the estimator itself may fail to satisfy asymptotic normality. In such cases, standard inference procedures—including CGM-type variance estimators as well as the approach proposed in this paper—become invalid. In this context, [31] proposes bootstrap-based inference methods that are uniformly valid, but they require the choice of tuning parameters and is typically conservative. In the two-way clustering setting, [32] develop a simple analytical inference that remains valid under non-Gaussian degeneracy, while avoiding the need for tuning parameters. Complementarily, [33] propose bootstrap procedures that are adaptive across a range of non-degenerate and (Gaussian) degenerate cases.

Extending these approaches to our setting is substantially more involved due to the presence of two-step estimators and machine learning-based first stages. In particular, degeneracy changes the effective stochastic order of the leading term, thereby tightening the rate requirements on the first-stage estimators to ensure that their estimation error is asymptotically negligible. Moreover, incorporating the strategy of [32] is highly non-trivial in our framework, as it relies on conditioning arguments that, in the presence of multiway dependence and generated regressors, further complicate the first-stage convergence requirements. Addressing these challenges would require new techniques, and we leave them for future research.”

4 Maximum Inequalities under separate exchangeability↩︎

In this section, we establish inequalities that control the \(q\)-th moment of the supremum of the empirical process, \(\mathbb{E}\bigl[\|\mathbb{G}_{n}\|_\mathcal{F}^q \bigr] ,\) for some \(q \in [1,\infty)\) for SE arrays. Throughout this section, assume without loss of generality that \(\mathbb{E}[f(X_{\boldsymbol{1}})] = 0\) for all \(f \in \mathcal{F}\). Before presenting the main results, let us first introduce the Hoeffding-type decomposition from [12]. For any \(\boldsymbol{i}\in [\boldsymbol{N}]\), define \[\begin{align} (P_{\boldsymbol{e}}f)\Bigl(\{U_{\boldsymbol{i}\odot \boldsymbol{e}'}\}_{\boldsymbol{e}'\le \boldsymbol{e}}\Bigr) &=\mathbb{E}\Bigl[f(X_{\boldsymbol{i}})\,\Big|\,\{U_{\boldsymbol{i}\odot \boldsymbol{e}'}\}_{\boldsymbol{e}'\le \boldsymbol{e}}\Bigr]. \end{align}\] We then define recursively for \(k=1,2,\dots, K\) that \[\begin{align} (\pi_{\boldsymbol{e}_k}f)(U_{\boldsymbol{i}\odot \boldsymbol{e}_k}) &=(P_{\boldsymbol{e}_k}f)(U_{\boldsymbol{i}\odot \boldsymbol{e}_k}), \end{align}\] and for \(\boldsymbol{e}\in \bigcup_{k=2}^K\mathcal{E}_k\) set \[\begin{align} (\pi_{\boldsymbol{e}}f)\Bigl(\{U_{\boldsymbol{i}\odot \boldsymbol{e}'}\}_{\boldsymbol{e}'\le \boldsymbol{e}}\Bigr) =\;& (P_{\boldsymbol{e}}f)\Bigl(\{U_{\boldsymbol{i}\odot \boldsymbol{e}'}\}_{\boldsymbol{e}'\le \boldsymbol{e}}\Bigr) \nonumber\\ &\quad-\sum_{\substack{\boldsymbol{e}'\le \boldsymbol{e}\\ \boldsymbol{e}'\ne \boldsymbol{e}}} (\pi_{\boldsymbol{e}'}f)\Bigl(\{U_{\boldsymbol{i}\odot \boldsymbol{e}''}\}_{\boldsymbol{e}''\le \boldsymbol{e}'}\Bigr). \end{align}\] Note that by the AHK representation 1 , for a fixed \(\boldsymbol{e}\) the distributions of \[(P_{\boldsymbol{e}}f)\Bigl(\{U_{\boldsymbol{i}\odot \boldsymbol{e}'}\}_{\boldsymbol{e}'\le \boldsymbol{e}}\Bigr) \quad\text{and}\quad (\pi_{\boldsymbol{e}}f)\Bigl(\{U_{\boldsymbol{i}\odot \boldsymbol{e}'}\}_{\boldsymbol{e}'\le \boldsymbol{e}}\Bigr)\] do not depend on the index \(\boldsymbol{i}\). Hence, we shall write \(P_{\boldsymbol{e}}f\) and \(\pi_{\boldsymbol{e}}f\) for a generic \(\boldsymbol{i}\).

Now, fix any \(1\le k\le K\) and let \(\boldsymbol{e}\in \mathcal{E}_k\). Then, by Lemma 1 in [12], for any \(\ell\in \mathrm{supp}(\boldsymbol{e})\) the random variable \((\pi_{\boldsymbol{e}}f)\Bigl(\{U_{\boldsymbol{i}\odot \boldsymbol{e}'}\}_{\boldsymbol{e}'\le \boldsymbol{e}}\Bigr)\) has mean zero conditionally on \(\{U_{\boldsymbol{i}\odot \boldsymbol{e}'}\}_{\boldsymbol{e}'\le \boldsymbol{e}-\boldsymbol{e}_{\ell}}\). In addition, define \(I_{\boldsymbol{N},\boldsymbol{e}} = \{\boldsymbol{i}\odot \boldsymbol{e} : \boldsymbol{i}\in [\boldsymbol{N}]\}.\) Then, we have \(\bigl|I_{\boldsymbol{N},\boldsymbol{e}}\bigr| = \prod_{k'\in\mathrm{supp}(\boldsymbol{e})} N_{k'}.\) Accordingly, define \[\begin{align} H_{\boldsymbol{N}}^{\boldsymbol{e}}(f) =\frac{1}{\bigl|I_{\boldsymbol{N},\boldsymbol{e}}\bigr|} \sum_{\boldsymbol{i}\in I_{\boldsymbol{N},\boldsymbol{e}}} (\pi_{\boldsymbol{e}}f)\Bigl(\{U_{\boldsymbol{i}\odot \boldsymbol{e}'}\}_{\boldsymbol{e}'\le \boldsymbol{e}}\Bigr). \end{align}\] We now obtain the Hoeffding-type decomposition \[\begin{align} \mathbb{E}_N f = \sum_{k=1}^K \sum_{\boldsymbol{e}\in\mathcal{E}_k} H_{\boldsymbol{N}}^{\boldsymbol{e}}(f). \label{hoeffding} \end{align}\tag{13}\] To bound \(\mathbb{E}\bigl[ \|\mathbb{G}_{n}(f)\|_{\mathcal{F}} \bigr]\), it thus suffices to control each individual term \(\mathbb{E}\bigl[\|H_{\boldsymbol{N}}^{\boldsymbol{e}}(f)\|_{\mathcal{F}}\bigr]\) separately.

Finally, fix any \(1\le k\le K\) and \(\boldsymbol{e}\in \mathcal{E}_k\). Define the uniform entropy integral by \[\begin{align} J_{\boldsymbol{e}}(\delta) = J_{\boldsymbol{e}}(\delta,\mathcal{F},F) := \int_{0}^{\delta} \sup_{Q} \Biggl\{ 1+ \log N\Bigl(P_{\boldsymbol{e}}\mathcal{F},\|\cdot\|_{Q,2},\tau\|P_{\boldsymbol{e}}F\|_{Q,2}\Bigr) \Biggr\}^{k/2} d\tau, \end{align}\] where \(P_{\boldsymbol{e}}\mathcal{F}:= \{P_{\boldsymbol{e}}f : f\in \mathcal{F}\},\) and the supremum is taken over all finite discrete distributions \(Q\).

The following result is a general global maximal inequality for SE empirical processes with an arbitrary index order \(K\) and for a general order of moment \(q\in[1,\infty)\). Its proof follows the arguments in the proof of Corollary B.1 in [12] with some modifications to account for a more general class of functions.

A proof can be found in Section 6.4 in the appendix.

Although the global maximal inequality works for general \(q\), in the case that the supremum of the first absolute moment is concerned, local maximal inequalities usually provides shaper bounds. The following is a novel local maximal inequality for SE empirical processes.

A proof can be found in Section 6.5 in the appendix.

In practice, bounding the uniform entropy integrals appearing on the right-hand side of maximal inequalities can be involved. Fortunately, many function classes arising in econometrics, machine learning, and statistics can be shown to be of VC-type, in the sense of 6 . Under this assumption, the entropy terms entering the maximal inequalities admit substantially simpler bounds. We therefore derive a local maximal inequality under the VC-type condition, which is the version utilised in the proof of Theorem 1.

A proof is provided in Section 6.6 in the appendix.

5 Conclusion↩︎

This paper develops a cross-fitting-free asymptotic theory for two-step debiased GMM estimators under multiway clustered dependence and shows that valid inference in such settings hinges on new empirical process techniques. By combining orthogonal moment conditions with a localisation-based argument, we demonstrate that the impact of high-dimensional or nonparametric nuisance estimation can be controlled without sample splitting, even when the effective sample size is determined by the number of independent cluster units. The resulting estimators are shown to be asymptotically linear and normal under separately exchangeable sampling, providing a practical inference framework for empirically relevant clustered environments. A central contribution is the derivation of new global and local maximal inequalities for possibly uncountable, pointwise measurable function classes under multiway dependence, which fill a gap in the existing theory and may be useful beyond the DML context. Future work may further explore stability-type conditions under multiway clustering and extend these tools to other two-step and high-dimensional problems with complex dependence structures.

Appendix↩︎

Throughout the appendix, for \(0 < \beta < \infty\), let \(\psi_{\beta}\) be the function on \([0,\infty)\) defined by \(\psi_{\beta} (x) = e^{x^{\beta}}-1\). Let \(\| \cdot \|_{\psi_{\beta}}\) denote the associated Orlicz norm, i.e., \(\| \xi \|_{\psi_\beta}=\inf \{ C>0: \mathbb{E}[ \psi_{\beta}( | \xi | /C)] \leq 1\}\) for a real-valued random variable \(\xi\).

6 Proofs of the main results↩︎

6.1 Proof of Theorem 1↩︎

Step 1. We consider the consistency of the Jacobian \(\widehat J_N(\widehat\theta)\). We decompose it as \(\Vert\widehat J_N(\widehat\theta) - J_0 \Vert \leq \Vert \widehat J_N(\widehat\theta) - \widehat J_N(\theta_0)\Vert + \Vert \widehat J_N(\theta_0) - J_0\Vert\). Consider \(\Vert \widehat J_N(\widehat\theta) - \widehat J_N(\theta_0)\Vert\): \[\begin{align} \Vert \widehat J_N(\widehat\theta) - \widehat J_N(\theta_0)\Vert =& \left\Vert \frac{1}{N} \sum_{\boldsymbol{i} \in [\boldsymbol{N}]} \partial_\theta\psi(X_{\boldsymbol{i}},\widehat\theta,\widehat\eta) - \partial_\theta\psi(X_{\boldsymbol{i}},\theta_0,\widehat\eta) \right\Vert \\ \leq & \frac{1}{N}\sum_{\boldsymbol{i} \in [\boldsymbol{N}]} \left\Vert \partial_\theta\psi(X_{\boldsymbol{i}},\widehat\theta,\widehat\eta) - \partial_\theta\psi(X_{\boldsymbol{i}},\theta_0,\widehat\eta) \right\Vert \\ \leq &\frac{1}{N}\sum_{\boldsymbol{i} \in [\boldsymbol{N}]} B(X_{\boldsymbol{i}},\widehat\eta) \Vert \widehat\theta-\theta_0\Vert^\alpha =O_p(1)o_P(1). \end{align}\] where the last line follows from Assumption 1(iii), \(\widehat\theta \overset{p}{\to}\theta\), and Lemma 6.

For \(\Vert \widehat J_N(\theta_0) - J_0\Vert\), we again apply Lemma 6 under the continuity of the score and \(\mathbb{E} [\sup_{\eta\in\Gamma_n(\eta_0)}\partial_\theta\psi(X,\theta_0,\eta)]<\infty\) ensured by Assumptions 1(ii) and (iii). So, we obtain \(\Vert\widehat J_N(\widehat\theta) - J_0 \Vert=o_P(1)\).

Step 2. Under the differentiability of the score with respect to \(\theta\), we can apply a mean-value expansion to the summand of \(\widehat\psi_N(\widehat\theta)\): \[\begin{align} \psi(X,\widehat\theta,\widehat\eta) = & \psi(X,\widehat\theta,\widehat\eta)- \psi(X,\theta_0,\widehat{\eta}) + \psi(X,\theta_0,\widehat{\eta}) - \psi(X,\theta_0,{\eta}_0) + \psi(X,\theta_0,{\eta}_0) \nonumber \\ =& \partial_\theta \widehat\psi(\theta)\big|_{\theta = \widetilde{\theta}}(\widehat\theta - \theta_0) + \psi(X,\theta_0,\widehat{\eta}) - \psi(X,\theta_0,{\eta}_0) + \psi(X,\theta_0,{\eta}_0) \label{eq:mv95expansion} \end{align}\tag{14}\] where \(\widetilde{\theta}\) is such that each of its coordinates lies between the corresponding element of \(\widehat\theta\) and \(\theta_0\). Then, plugging (14 ) to (5 ) gives \[\begin{align} \sqrt{n}(\widehat\theta - \theta_0) = & \left( \widehat J_N(\widehat\theta)' \widehat\Upsilon \widehat J_N(\widetilde{\theta})\right)^{-1} \widehat J_N(\widehat\theta)' \widehat\Upsilon \sqrt{n} \mathbb{E}_N [\psi(X,\theta_0,{\eta}_0)] \\ &+\left( \widehat J_N(\widehat\theta)' \widehat\Upsilon \widehat J_N(\widetilde{\theta})\right)^{-1} \widehat J_N(\widehat\theta)' \widehat\Upsilon \mathbb{G}_{n}(f(\widehat\eta))\\ &+\left( \widehat J_N(\widehat\theta)' \widehat\Upsilon \widehat J_N(\widetilde{\theta})\right)^{-1} \widehat J_N(\widehat\theta)' \widehat\Upsilon \sqrt{n} \mathbb{E}[ f(\eta)]_{\eta = \widehat\eta}, \end{align}\] where, recall that, \(\mathbb{G}_{n}\left(f(\eta)\right): =\sqrt{n}\left( \mathbb{E}_N [ f(\eta)] - \mathbb{E}[ f(\eta)]\right)\) and \(f(\eta) = \psi(X,\theta_0,{\eta}) - \psi(X,\theta_0,{\eta}_0)\). It can be shown similarly that \(\Vert \widehat J_N(\widetilde{\theta}) - J_0 \Vert = o_P(1)\) as above. Given the consistency of \(\widehat\Upsilon\) for some positive-definite \(\Upsilon\) and that \({\rm rank}(J_0)\) = d, \(J_0'\Upsilon J_0\) is nonsingular. It is then left to show a CLT for \(\mathbb{E}_N [\psi(X,\theta_0,{\eta}_0)]\) and to bound \(\mathbb{G}_{n}(f(\widehat\eta))+ \sqrt{n} \mathbb{E}[ f(\eta)]_{\eta = \widehat\eta}\).

Step 2-1. We first bound \(\mathbb{E}[ f(\eta)]_{\eta = \widehat\eta}\) using the orthogonality condition. Due to Gateaux differentiability by Assumption 1(ii), the fundamental theorem of calculus, and the Neyman orthogonality condition (3 ), \[\begin{align} &\left\Vert \mathbb{E}[ f(\eta)]_{\eta = \widehat\eta} \right\Vert = \left\Vert \int_0^1 \partial_{\tau} \mathbb{E}\left[\psi(X,\theta_0,\eta_0 + \tau(\widehat\eta-\eta_0))\right] d\tau\right\Vert \\ =&\left\Vert \int_0^1 \left(\partial_{\tau}\mathbb{E}\left[\psi(X,\theta_0,\eta_0 + \tau(\widehat\eta-\eta_0))\right] - \partial_{\tau}\mathbb{E}\left[\psi(X,\theta_0,\eta_0 + \tau(\widehat\eta-\eta_0))\right]_{\tau=0} \right)d\tau\right\Vert \\ \leq& \int_0^1 \left\Vert\partial_{\eta}\mathbb{E}\left[\psi(X,\theta_0,\eta_0 + \tau(\widehat\eta-\eta_0))\right] - \partial_{\eta}\mathbb{E}\left[\psi(X,\theta_0,\eta_0))\right] d\tau\right\Vert \left\Vert \widehat\eta-\eta_0\right\Vert_{P,2} \end{align}\] Furthermore, by the Lipschitz continuity of the Gateaux derivative, there exists some constant \(C_4<\infty\) such that \(\left\Vert\partial_{\eta}\mathbb{E}\left[\psi(X,\theta_0,\eta_0 + \tau(\widehat\eta-\eta_0))\right] - \partial_{\eta}\mathbb{E}\left[\psi(X,\theta_0,\eta_0))\right] d\tau\right\Vert \leq C_4 \tau\Vert\widehat\eta-\eta_0 \Vert_{P,2}\). Therefore, under Assumption 1(iv), we have \[\begin{align} \Vert \sqrt{n} \mathbb{E}[ f(\eta)]_{\eta=\widehat\eta}\Vert = \frac{\sqrt{n}C_4}{2} \Vert\widehat\eta-\eta_0 \Vert_{P,2}^2 = o_P(1). \end{align}\]

Step 2-2. Next, we will bound \(\mathbb{G}_{n}(f(\widehat\eta))\) through the maximum inequality for multiway clustering data. Recall that following Assumption 1, there is a local neighborhood \(\Gamma_n(\eta_0)= \{\eta: \Vert \eta-\eta_0\Vert_{P,2} \le C_1n^{-1/4}\}\). Then we have \(\widehat\eta\in \Gamma_n(\eta_0)\) with probability converging to one. We also define a class of functions for the centered \(f(\eta)\): \[\begin{align} \mathcal{F}_{n}: = \{ f(\eta) - \mathbb{E}[ f(\eta)]:\eta\in \Gamma_n(\eta_0)\}. \end{align}\]

Recall that \({H}_{\boldsymbol{N}}^{\boldsymbol{e}}(f) =\bigl|I_{\boldsymbol{N},\boldsymbol{e}}\bigr|^{-1} \sum_{\boldsymbol{i}\in I_{\boldsymbol{N},\boldsymbol{e}}} \pi_{\boldsymbol{e}}f\Bigl(\{U_{\boldsymbol{i}\odot \boldsymbol{e}'}\}_{\boldsymbol{e}'\le \boldsymbol{e}}\Bigr)\) where \(I_{\boldsymbol{N},\boldsymbol{e}} = \{\boldsymbol{i}\odot \boldsymbol{e} : \boldsymbol{i}\in [\boldsymbol{N}]\}\) and \(\bigl|I_{\boldsymbol{N},\boldsymbol{e}}\bigr| = \prod_{k'\in\mathrm{supp}(\boldsymbol{e})} N_{k'}.\). By Hoeffding decomposition given in Section 4 and the triangle inequality, we have with probability approaching one \[\begin{align} \Vert \mathbb{G}_{n}\left(f(\widehat\eta)\right)\Vert\leq \sup_{\eta\in \Gamma_n} \left\Vert\sqrt{n} \mathbb{E}_N \left[ f(\eta) - \mathbb{E}[ f(\eta)]\right]\right\Vert \leq \sqrt{n} \sum_{k=1}^K \sum_{\boldsymbol{e}\in\mathcal{E}_k}\sup_{f\in \mathcal{F}_{n}}\left\Vert {H}_{\boldsymbol{N}}^{\boldsymbol{e}}(f)\right\Vert \end{align}\]

For each \(\boldsymbol{e}\in\mathcal{E}_k\) and each \(k=1,...,K\), we can apply Theorem 5: \[\sqrt{n}\sup_{f\in \mathcal{F}_{n}}\left\Vert {H}_{\boldsymbol{N}}^{\boldsymbol{e}}(f)\right\Vert\lesssim \sqrt{\frac{n}{|I_{\boldsymbol{N,e}}|}} \left(J_{\boldsymbol{e}}(\delta_{\boldsymbol{e}}) \Vert P_{\boldsymbol{e}} F_n \Vert_{P,2} + \frac{J^2_{\boldsymbol{e}}(\delta_{\boldsymbol{e}})\Vert M_{\boldsymbol{e}}\Vert_{P,2}}{\sqrt{n}\delta^2_{\boldsymbol{e}}}\right).\] Combining with the consistency of \(\widehat{\Upsilon}_N\) and \(\widehat J_N(\widehat\theta)\overset{p}{\to}\) shown above, we obtain the following: \[\begin{align} \sqrt{n}(\widehat{\theta}-\theta_0) = &\left( J_0' \Upsilon J_0\right)^{-1} J_0' \Upsilon \sqrt{n} \mathbb{E}_N [\psi(X,\theta_0,{\eta}_0)] + \sum_{k=1}^K\sum_{\boldsymbol{e}\in \mathcal{E}_k} O_P(\rho_{n,\boldsymbol{e}}) +o_P(1) \end{align}\] where \(\rho_{n,\boldsymbol{e}}=\sqrt{\frac{n}{|I_{\boldsymbol{N},e}|}} \left(J_{\boldsymbol{e}}(\delta_{\boldsymbol{e}}) \Vert P_e F_n \Vert_{P,2} + \frac{J^2_{\boldsymbol{e}}(\delta_{\boldsymbol{e}})\Vert M_{\boldsymbol{e}}\Vert_{P,2}}{\sqrt{n}\delta^2_{\boldsymbol{e}}}\right)\).

Step 3. For the first statement, under Assumption 1(ii), we have \[\begin{align} \sup_{\eta\in \Gamma_n} \mathbb{E}\|f(\eta)\|^2 = \sup_{\eta\in \Gamma_n} \mathbb{E} \|\psi(X,\theta_0,{\eta}) - \psi(X,\theta_0,{\eta}_0)\|^2 \leq \sup_{\eta\in \Gamma_n}C_2\Vert \eta - \eta_0\Vert^2_{P,2}= O(n^{-1/2}) \end{align}\] It also follows that \(\sup_{\eta\in \Gamma_n} \mathbb{E}\|f(\eta) - \mathbb{E}[ f(\eta)]\|^2 =O(n^{-1/2})\). We observe that \({H}_{\boldsymbol{N}}^{\boldsymbol{e}}(f)\) is a linear combination of \(P\boldsymbol{e}(f)\). By Jensen’s inequality and law of iterated expectation, we have \[\begin{align} \sup_{f\in\mathcal{F}_{n}} \mathbb{E}\|P\boldsymbol{e}(f)\|^2 \lesssim \sup_{\eta\in \Gamma_n} \mathbb{E}\|f(\eta) - \mathbb{E}[ f(\eta)]\|^2 =O(n^{-1/2}). \end{align}\]

Therefore, we can take \(\sigma_{\boldsymbol{e}}=\sigma_{{\boldsymbol{e}},n}\) as a sequence of positive numbers with \(\sigma_{\boldsymbol{e}} = O(n^{-1/4})\). If \(\mathcal{F}_{n}\) is a VC-type class with characteristics \(A_n\ge (e^{2(K-1)/16}\vee e)\) and \(v\ge1\), then applying Corollary \(1\) gives, for all \(\boldsymbol{e}\in \mathcal{E}_k\) and \(k=1,...,K\) \[\begin{align} & \sqrt{n} E\left[\sup_{f\in\mathcal{F}_{n}}\|{H}_{\boldsymbol{N}}^{\boldsymbol{e}}(f)\|\right] \\ =& \sqrt{\frac{n}{|I_{\boldsymbol{N},\boldsymbol{e}}|}} |I_{\boldsymbol{N},\boldsymbol{e}}|^{1/2} \mathbb{E}\left[\sup_{f\in\mathcal{F}_{n}}\|{H}_{\boldsymbol{N}}^{\boldsymbol{e}}(f)\|\right] \\ \lesssim& \frac{n^{1/4}}{|I_{\boldsymbol{N},\boldsymbol{e}}|^{1/2}} \{v_n\log(A_n\vee \overline{N} )\}^{k/2} + \frac{\|M_{\boldsymbol{e}}\|_{P,2}}{|I_{\boldsymbol{N},\boldsymbol{e}}|^{1/2}} \{v_n\log(A_n\vee \overline{N} )\}^{k} . \end{align}\] where \(M_{\boldsymbol{e}}=\max_{t\in I_{\boldsymbol{N},\boldsymbol{e}}}(P_{\boldsymbol{e}}F_n)\left(\{U_{(t,...,t)\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}\right)\), and, for some \(q>2\), \[\begin{align} &\|M_{\boldsymbol{e}}\|_{P,2} \leq \|M_{\boldsymbol{e}}\|_{P,q} \leq \left\| \sum_{t\in I_{\boldsymbol{N},\boldsymbol{e}}} (P_{\boldsymbol{e}}F_n)\left(\{U_{(t,...,t)\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}\right) \right\|_{P,q} \\ \leq &\left( \sum_{t\in I_{\boldsymbol{N},\boldsymbol{e}}}\mathbb{E} \left\| (P_{\boldsymbol{e}}F_n)\left(\{U_{(t,...,t)\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}\right) \right\|^q \right)^{1/q}\leq \left|I_{\boldsymbol{N},\boldsymbol{e}}\right|^{1/q} \left\| F_n \right\|_{P,q}. \end{align}\] Then, by Markov’s inequality, \[\begin{align} \sqrt{n}\sup_{f\in \mathcal{F}_{n}}\left\Vert {H}_{\boldsymbol{N}}^{\boldsymbol{e}}(f)\right\Vert \lesssim& \sqrt{ \frac{\{v_n\log(A_n\vee \overline{N} )\}^{k}}{|I_{\boldsymbol{N},\boldsymbol{e}}|}} + \frac{\left\| F_n \right\|_{P,q} \{v_n\log(A_n\vee \overline{N} )\}^{k}}{|I_{\boldsymbol{N},\boldsymbol{e}}|^{1/2-1/q}} \\ \leq &\left(\frac{v_n\log(A_n\vee \overline{N} )}{n^{1-1/2k}}\right)^{k/2} + \left(\frac{\left\| F_n \right\|_{P,q} \{v_n\log(A_n\vee \overline{N} )\}}{n^{1/2-1/q}}\right)^{k}, \end{align}\] which proves the first statement.

For the second statement, we apply Lemma 7 to obtain an independent linear representation, \[\begin{align} &\sqrt{n}\mathbb{E}_N[\psi(X_{\boldsymbol{i}},\theta_0,\eta_0)] = \sum_{\boldsymbol{i} \in I_1} \frac{\sqrt{n}}{N_{k(\boldsymbol{i})}} \mathbb{E}\left[\psi(X_{\boldsymbol{i}},\theta_0,\eta_0) | U_{\boldsymbol{i}}\right] + O_P(n^{-1/2}) \nonumber \\ =& \frac{\sqrt{n}}{N_1}\sum_{i_1=1}^{N_1} \mathbb{E}[\psi(X_{\boldsymbol{i}},\theta_0,\eta_0)|U_{1,0,...,0}] + ...+ \frac{\sqrt{n}}{N_K}\sum_{i_K=1}^{N_K} \mathbb{E}[\psi(X_{\boldsymbol{i}},\theta_0,\eta_0)|U_{0,0,...,1}] +O_P(n^{-1/2}) \label{eq:hajek} \end{align}\tag{15}\] and the variance of \(\sqrt{n}\mathbb{E}_N[\psi(X_{\boldsymbol{i}},\theta_0,\eta_0)]\), \[\begin{align} & {\rm Var} \left( \sqrt{n}\mathbb{E}_N[\psi(X_{\boldsymbol{i}}, \theta_0,\eta_0)] \right) =\sum_{\boldsymbol{e}\in\mathcal{E}_1} \mu_{k(\boldsymbol{e})} {\rm Cov}(\psi(X_{\boldsymbol{1}},\theta_0,\eta_0),\psi(X_{\boldsymbol{2}-\boldsymbol{e}},\theta_0,\eta_0)) + O(n^{-1}) \nonumber \\ =& \mu_1(1+o(1)) {\rm Var}( \mathbb{E}[\psi(X,\theta_0,\eta_0)|U_{1,0,...,0}]) + ... +\mu_K(1+o(1)) {\rm Var}( \mathbb{E}[\psi(X,\theta_0,\eta_0)|U_{0,0,...,1}]) \\&+ O(n^{-1}), \end{align}\] where \(k(\boldsymbol{e})\), denote the coordinate in which \(\boldsymbol{e}\in \mathcal{E}_1\) is non-zero.

Note that each term \(\mathbb{E}\left[\psi(X_{\boldsymbol{i}},\theta_0,\eta_0) | U_{\boldsymbol{i}}\right]\) is independently and identically distributed. Given the finite second moment of \(\psi(X_{\boldsymbol{i}},\theta_0,\eta_0)\) in Assumption 1(v) (and so ) and that \(\psi(X_{\boldsymbol{i}},\theta_0,\eta_0)\) is mean-zero, we can apply Lindeberg–Lévy CLT as well as Cramer-Wold device to each of the sums in (15 ) and obtain \[\begin{align} \sqrt{n}\mathbb{E}_N[\psi(X_{\boldsymbol{i}},\theta_0,\eta_0)] \overset{d}{\to} N\left(0, \lim_{n\to \infty}{\rm Var} \left( \sqrt{n}\mathbb{E}_N[\psi(X_{\boldsymbol{i}}, \theta_0,\eta_0)] \right) \right) \end{align}\] Given the rate condition of 8 , the second statement follows. 0◻

6.2 Proof of Theorem 2↩︎

Case (1). For any index set \(I\subset{1,...,p}\) such that \(|I| =s\), we can define a subclass of \(\mathcal{G}_n\) as \(\mathcal{G}_I^* = \{x\mapsto x^T_I\beta_I:\beta_I\in \mathbb{R}^{|I|},\Vert\beta_I \Vert_2<\infty \}\). Let \(G\) be an envelope of \(\mathcal{G}_n\). Then, for any finite discrete measure \(Q\) and for some \(C<\infty\), \[\begin{align} N(\mathcal{G}_I,\Vert \cdot \Vert_{Q,2},\varepsilon,\Vert G\Vert_{Q,2}) \leq \left(\frac{C}{\varepsilon}\right)^s. \end{align}\] We note that \(\mathcal{G}_I^* = \bigcup_{|I|=s}\mathcal{G}_I\). Therefore, we can bound the covering number as follows: \[\begin{align} N(\mathcal{G}_n,\Vert \cdot \Vert_{Q,2},\varepsilon,\Vert G\Vert_{Q,2}) \leq& \binom{p}{s} \left(\frac{C}{\varepsilon}\right)^s \leq \left(\frac{ep}{s}\right)^s\left(\frac{C}{\varepsilon}\right)^s = \left(\frac{C \;ep}{s\varepsilon}\right)^s. \end{align}\] Therefore, we can choose \(v_n= s\) and \(A_n = \frac{C \;ep}{s}\vee (e^{2(K-1)/16}\vee e)\). Given this choice and that \({s\log(p/s)} \vee s\log(\overline{N})= o(n^{1/4})\), we have \(\frac{v_n \log(A_n\vee \overline{N})}{n^{1-1/2k}} = o(1)\) for all \(k\geq 1\) and \(\frac{v_n \log(A_n\vee \overline{N})}{n^{1/2-1/q}} = o(1)\) for some \(q\geq 4\). The statement of the case (1) follows by noting that the monotonic transformation preserves the VC-type characteristics.

Case (2). Fix any \(g\in \mathcal{G}_n\) and \(r\in\mathbb{R}\). We can build a decision tree by setting the splitting rule \(t(x) = 1\{g(x)>r\}\), and the corresponding class can be written as \[\begin{align} \mathcal{H}_N = \{(x,r)\mapsto 1\{g(x)>r\}: g\in \mathcal{G}_n, \;r\in\mathbb{R}\} \end{align}\] We note that the VC dimension of the subgraph of \(\mathcal{G}_n\) is equal to \({\rm VC}(\mathcal{H}_N)\), and the decision tree associated with \(\mathcal{H}_N\) has \(2L\) partitions.

Due to [35], the class of classification functions induced by a binary decision tree with partition \(\{R_l\}_{l=1}^{2L}\) defined in the statement has VC dimension of order \(O(L\log(2Lp))\) where \(p\) is the dimension of the real-valued features. Therefore, \(\mathcal{G}_n\) is a VC-subgraph class with the index \(V =CL\log(2Lp)\) for some constant \(C\in (0,\infty)\). Let \(G\) be the envelope of \(\mathcal{G}_n\). Since \(\mu_l<\infty\) for all \(l\), we can take \(G\) as some large enough constant. By Theorem 2.6.7 in [36], \[\begin{align} & N(\mathcal{G}_n,\Vert \cdot \Vert_{Q,2},\varepsilon,\Vert G\Vert_{Q,2}) \leq K V(16e)^V (1/\varepsilon)^{2(V-1)}\\ \lesssim & CK L\log(2Lp) (16e/\varepsilon)^{2CL\log(2Lp)} \\ =& K \left([CL\log(2Lp)]^{1/2CL\log(2Lp)}16e/\varepsilon\right)^{2CL\log(2Lp)} \end{align}\] for a universal constant \(K\). Therefore, we can take \(v_n = 2CL\log(2Lp)\) and \(A_n\) as some large enough constant because \([CL\log(2Lp)]^{1/2CL\log(2Lp)}\) is bounded for fixed \(L\) and converges to some constant as \(L\) diverges. The rate condition follows immediately.

Case (3). Let \(\mathcal{G}_n^{+}\) denote the class of functions induced by the neural network defined in the statement but with \(L-1\) hidden layers and \(W\) parameters. Let \({\rm sgn}(\mathcal{G}_{N}) : = \{{\rm sgn}(g): g\in \mathcal{G}_{N}\}\) and \({\rm sub}(\mathcal{G}_{N}) = \{(x,y)\in \mathcal{X}\times\mathbb{R}: g(x)>y, g\in \mathcal{G}_n\}\). By Theorem 14.1 of [37], \({\rm VC}({\rm sub}(\mathcal{G}_{N}))\leq{\rm VC}({\rm sgn}(\mathcal{G}_{N}^{+})\). By Theorem 7 of [38], \({\rm VC}({\rm sgn}(\mathcal{G}_{N}^{+})) = O(LW\log(pU))\). Then, by Theorem 2.6.7 of [36] and a similar calculation of Case (2), we can choose \(A_n= C\vee(e^{2(K-1)/16}\vee e)\) for some large enough constant \(C\) and \(v_n = 2LW\log(pU)\). The rate condition then follows.

0◻

6.3 Proof of Theorem 3↩︎

6.4 Proof of Theorem 4↩︎

By symmetrisation inequality for SE processes (Lemma B.1 in [12]; note that it is dimension free), for independent Rademacher r.v.’s \((\varepsilon_{1,i_1})\),...,\((\varepsilon_{k,i_k})\) that are independent of \((X_{\boldsymbol{i}})_{{\boldsymbol{i}}\in\mathbb{N}^K}\), one has \[\begin{align} |I_{{\boldsymbol{N}},{\boldsymbol{e}}}|^{1/2} \left(\mathbb{E}[\|H_{\boldsymbol{N}}^{\boldsymbol{e}}(f)\|_\mathcal{F}^q]\right)^{1/q} =&\left (\mathbb{E} \left [ \left \|\frac{1}{\sqrt{|I_{{\boldsymbol{N}},{\boldsymbol{e}}}|}} \sum_{{\boldsymbol{i}}\in I_{{\boldsymbol{N}},{\boldsymbol{e}}}} (\pi_{ {\boldsymbol{e}}}f )(\{U_{{\boldsymbol{i}}\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}) \right \|_{\mathcal{F}}^{q} \right ] \right)^{1/q}\\ \lesssim& \left (\mathbb{E} \left [ \left \|\frac{1}{\sqrt{|I_{{\boldsymbol{N}},{\boldsymbol{e}}}|}} \sum_{{\boldsymbol{i}}\in I_{{\boldsymbol{N}},{\boldsymbol{e}}}}\varepsilon_{1,i_1}...\varepsilon_{k,i_k}\cdot(\pi_{ {\boldsymbol{e}}}f )(\{U_{{\boldsymbol{i}}\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}) \right \|_{\mathcal{F}}^{q} \right ] \right)^{1/q}. \end{align}\] By convexity of supremum and \(\cdot \mapsto (\cdot)^q\), Jensen’s inequality implies that the RHS above can be upperbounded up to a constant that depends only on \(q\), \(K\), and \(k\) by \[\begin{align} \left (\mathbb{E} \left [ \left \|\frac{1}{\sqrt{|I_{{\boldsymbol{N}},{\boldsymbol{e}}}|}} \sum_{{\boldsymbol{i}}\in I_{{\boldsymbol{N}},{\boldsymbol{e}}}}\varepsilon_{1,i_1}...\varepsilon_{k,i_k}\cdot(P_{\boldsymbol{e}}f )(\{U_{{\boldsymbol{i}}\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}) \right \|_{\mathcal{F}}^{q} \right ] \right)^{1/q}. \end{align}\] Denote \(\mathbb{P}_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}}=|I_{{\boldsymbol{N}},{\boldsymbol{e}}}|^{-1}\sum_{{\boldsymbol{i}}\in I_{{\boldsymbol{N}},{\boldsymbol{e}}}}\delta_{\{U_{{\boldsymbol{i}}\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}}\), the empirical measure on the support of \(\{U_{{\boldsymbol{i}}\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}\). Observe that conditionally on \(\{X_{\boldsymbol{i}}\}_{{\boldsymbol{i}}\in{\boldsymbol{N}}}\), the object \[\begin{align} \frac{1}{\sqrt{|I_{{\boldsymbol{N}},{\boldsymbol{e}}}|}}\sum_{{\boldsymbol{i}}\in I_{{\boldsymbol{N}},{\boldsymbol{e}}}}\varepsilon_{1,i_1}...\varepsilon_{k,i_k}(P_{{\boldsymbol{e}}}f )(\{U_{{\boldsymbol{i}}\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}) \end{align}\] is a homogeneous Rademacher chaos process of order \(k\). By Lemma 3, \(L^q\) norm is bounded from above by \(\psi_{2/k}\)-norm up to a constant depends only on \((q,k)\), and thus by applying Corollary 5,1.8 in [39], one has \[\begin{align} &\left (\mathbb{E} \left [ \left \|\frac{1}{\sqrt{|I_{{\boldsymbol{N}},{\boldsymbol{e}}}|}} \sum_{{\boldsymbol{i}}\in I_{{\boldsymbol{N}},{\boldsymbol{e}}}}\varepsilon_{1,i_1}...\varepsilon_{k,i_k} \cdot(P_{ {\boldsymbol{e}}}f )(\{U_{{\boldsymbol{i}}\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}) \right \|_{\mathcal{F}}^{q} \right ] \right)^{1/q}\\ \lesssim& \mathbb{E} \left [ \left\| \left \|\frac{1}{\sqrt{|I_{{\boldsymbol{N}},{\boldsymbol{e}}}|}} \sum_{{\boldsymbol{i}}\in I_{{\boldsymbol{N}},{\boldsymbol{e}}}}\varepsilon_{1,i_1}...\varepsilon_{k,i_k}\cdot(P_{ {\boldsymbol{e}}}f )(\{U_{{\boldsymbol{i}}\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}) \right \|_{\mathcal{F}}\right\|_{\psi_{2/k}|(X_{\boldsymbol{i}})_{{\boldsymbol{i}}\in [{\boldsymbol{N}}]}}\right ] \\ \lesssim& \mathbb{E} \left [ \int_0^{\sigma_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}}}\left[1+\log N\left(P_{ {\boldsymbol{e}}}\mathcal{F},\|\cdot\|_{\mathbb{P}_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}},2},\tau\right)\right]^{k/2}d \tau \right ], \end{align}\] where \(\sigma_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}}^2:=\sup_{f\in\mathcal{F}}\| P_{\boldsymbol{e}}f \|_{\mathbb{P}_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}},2}^2\). Using a change of variable and the definition of \(J_{\boldsymbol{e}}\), the above bound becomes \[\begin{align} &\mathbb{E} \left [ \int_0^{\sigma_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}}}\left[1+\log N\left(P_{ {\boldsymbol{e}}} \mathcal{F},\|\cdot\|_{\mathbb{P}_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}},2},\tau\right)\right]^{k/2}d \tau \right ] \\ =&\mathbb{E} \left [\|P_{ {\boldsymbol{e}}} F\|_{\mathbb{P}_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}},2} \int_0^{\sigma_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}}/\|P_{ {\boldsymbol{e}}} F\|_{\mathbb{P}_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}},2}}\left[1+\log N\left(P_{ {\boldsymbol{e}}} \mathcal{F},\|\cdot\|_{\mathbb{P}_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}},2},\tau\|P_{ {\boldsymbol{e}}} F\|_{\mathbb{P}_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}},2}\right)\right]^{k/2}d \tau \right ]\\ \le& \mathbb{E} \left [\|P_{ {\boldsymbol{e}}} F\|_{\mathbb{P}_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}},2} J_{\boldsymbol{e}}\left(\sigma_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}}/\|P_{ {\boldsymbol{e}}} F\|_{\mathbb{P}_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}},2}\right) \right ]\\ \le&J_{\boldsymbol{e}}\left(1\right) \|F\|_{P,q\vee 2} , \end{align}\] where the last inequality follows from Jensen’s inequality. 0◻

6.5 Proof of Theorem 5↩︎

We first state a crucial technical lemma which will be used in the following proof, a proof of this lemma is provided in the end of this section.

We now present the proof of Theorem 5. For an \({\boldsymbol{e}}\in \mathcal{E}_1\), the summands are i.i.d. and thus the desired result follows directly from Lemma 2. Therefore, we assume \(K\ge 2\) and \({\boldsymbol{e}}\in \mathcal{E}_k\) for a \(k\in\{2,...,K\}\). Assume without loss of generality that \({\boldsymbol{e}}\) consists of \(1\)’s in its first \(k\) elements and zero elsewhere. By applying the symmetrisation of Lemma B.1 in [12], one has, for independent Rademacher r.v.’s \((\varepsilon_{1,i_1})\),...,\((\varepsilon_{k,i_k})\) that are independent of \((X_{\boldsymbol{i}})_{{\boldsymbol{i}}\in\mathbb{N}^K}\), that \[\begin{align} |I_{{\boldsymbol{N}},{\boldsymbol{e}}}|^{1/2}\mathbb{E}[\|H_{\boldsymbol{N}}^{\boldsymbol{e}}(f)\|_\mathcal{F}] \lesssim& \mathbb{E}\left[\left|\frac{1}{\sqrt{|I_{{\boldsymbol{N}},{\boldsymbol{e}}}|}}\sum_{{\boldsymbol{i}}\in I_{{\boldsymbol{N}},{\boldsymbol{e}}}}\varepsilon_{1,i_1}...\varepsilon_{k,i_k}(\pi_{{\boldsymbol{e}}}f )(\{U_{{\boldsymbol{i}}\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}})\right\|_\mathcal{F}\right]. \end{align}\] Further, by convexity of supremum and Jensen’s inequality, the RHS above can be upper-bounded up to a constant that depends only on \(K\) and \(k\) by \[\begin{align} \mathbb{E}\left[\left|\frac{1}{\sqrt{|I_{{\boldsymbol{N}},{\boldsymbol{e}}}|}}\sum_{{\boldsymbol{i}}\in I_{{\boldsymbol{N}},{\boldsymbol{e}}}}\varepsilon_{1,i_1}...\varepsilon_{k,i_k}(P_{{\boldsymbol{e}}}f )(\{U_{{\boldsymbol{i}}\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}})\right\|_\mathcal{F}\right] \end{align}\]

Denote \(\mathbb{P}_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}}=|I_{{\boldsymbol{N}},{\boldsymbol{e}}}|^{-1}\sum_{{\boldsymbol{i}}\in I_{{\boldsymbol{N}},{\boldsymbol{e}}}}\delta_{\{U_{{\boldsymbol{i}}\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}}\), the empirical measure on the support of \(\{U_{{\boldsymbol{i}}\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}\). Observe that conditionally on \(\{X_{\boldsymbol{i}}\}_{{\boldsymbol{i}}\in{\boldsymbol{N}}}\), the object \[\begin{align} R_{\boldsymbol{N}}^{\boldsymbol{e}}(f)=\frac{1}{\sqrt{|I_{{\boldsymbol{N}},{\boldsymbol{e}}}|}}\sum_{{\boldsymbol{i}}\in I_{{\boldsymbol{N}},{\boldsymbol{e}}}}\varepsilon_{1,i_1}...\varepsilon_{k,i_k}(P_{{\boldsymbol{e}}}f )(\{U_{{\boldsymbol{i}}\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}) \end{align}\] is a homogeneous Rademacher chaos process of order \(k\). Further, following Corollary 3.2.6 in [39], for any \(f,f'\in \mathcal{F}\) \[\begin{align} \left\|R_{\boldsymbol{N}}^{\boldsymbol{e}}(f)-R_{\boldsymbol{N}}^{\boldsymbol{e}}(f')\right\|_{\psi_{2/k}|\{X_{\boldsymbol{i}}\}_{{\boldsymbol{i}}\in {\boldsymbol{N}}}}\lesssim \left\|R_{\boldsymbol{N}}^{\boldsymbol{e}}(f)-R_{\boldsymbol{N}}^{\boldsymbol{e}}(f')\right\|_{\mathbb{P}_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}},2}. \end{align}\] Hence the diameter of the function class \(\mathcal{F}\) in \(\|\cdot\|_{\psi_{2/k}|\{X_{\boldsymbol{i}}\}_{{\boldsymbol{i}}\in {\boldsymbol{N}}}}\)-norm is upperbounded by \(\sigma_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}}^2\) up to a constant, where \(\sigma_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}}^2:=\sup_{f\in\mathcal{F}}\| P_{{\boldsymbol{e}}}f \|_{\mathbb{P}_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}},2}^2\). By applying Fubini’s theorem, Corollary 5,1.8 in [39], and a change of variables, we have \[\begin{align} &\mathbb{E}\left[\left\|\frac{1}{\sqrt{|I_{{\boldsymbol{N}},{\boldsymbol{e}}}|}}\sum_{{\boldsymbol{i}}\in I_{{\boldsymbol{N}},{\boldsymbol{e}}}}\varepsilon_{1,i_1}...\varepsilon_{k,i_k}(P_{{\boldsymbol{e}}}f )(\{U_{{\boldsymbol{i}}\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}})\right\|_\mathcal{F}\right]\\ \lesssim& \mathbb{E}\left[\left\|\left\|\frac{1}{\sqrt{|I_{{\boldsymbol{N}},{\boldsymbol{e}}}|}}\sum_{{\boldsymbol{i}}\in I_{{\boldsymbol{N}},{\boldsymbol{e}}}}\varepsilon_{1,i_1}...\varepsilon_{k,i_k}(P_{{\boldsymbol{e}}}f )(\{U_{{\boldsymbol{i}}\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}})\right\|_\mathcal{F}\right\|_{\psi_{2/k}|\{X_{\boldsymbol{i}}\}_{{\boldsymbol{i}}\in{\boldsymbol{N}}}}\right]\\ \lesssim& \mathbb{E} \left [ \int_0^{\sigma_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}}}\left[1+\log N\left(P_{ {\boldsymbol{e}}}\mathcal{F},\|\cdot\|_{\mathbb{P}_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}},2},\tau\right)\right]^{k/2}d \tau \right ] \\ =&\mathbb{E} \left [ \|P_{\boldsymbol{e}}F\|_{\mathbb{P}_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}},2} \int_0^{\sigma_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}}/\|P_{\boldsymbol{e}}F\|_{\mathbb{P}_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}},2}}\left[1+\log N\left(P_{ {\boldsymbol{e}}}\mathcal{F},\|\cdot\|_{\mathbb{P}_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}},2},\tau\|P_{\boldsymbol{e}}F\|_{\mathbb{P}_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}},2}\right)\right]^{k/2}d \tau \right ]\\ \le& \mathbb{E} \left [\|P_{\boldsymbol{e}}F\|_{\mathbb{P}_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}},2} J_{\boldsymbol{e}}\left(\sigma_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}}/\|P_{\boldsymbol{e}}F\|_{\mathbb{P}_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}},2}\right) \right ]. \end{align}\] By Lemma 4, an application of Jensen’s inequality yields \[\begin{align} |I_{{\boldsymbol{N}},{\boldsymbol{e}}}|^{1/2}\mathbb{E}[\|H_{{\boldsymbol{N}}}^{\boldsymbol{e}}(f)\|_\mathcal{F}]\lesssim& \|P_{\boldsymbol{e}}F\|_{P,2} J_{\boldsymbol{e}}\left(z\right), \label{eq:local95maximal95ineq95prelimary95bound} \end{align}\tag{16}\] where \(z:=\sqrt{\mathbb{E}[\sigma_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}}^2]/\|P_{\boldsymbol{e}}F\|_{P,2}^2}\).

We now bound \[\begin{align} \mathbb{E}[\sigma_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}}^2]=&\mathbb{E}\left[\left\|\frac{1}{|I_{{\boldsymbol{N}},{\boldsymbol{e}}}|}\sum_{{\boldsymbol{i}}\in I_{{\boldsymbol{N}},{\boldsymbol{e}}}}(P_{\boldsymbol{e}}f )^2\left(\{U_{{\boldsymbol{i}}\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}\right)\right\|_\mathcal{F}\right]. \end{align}\] We aim to apply the Hoffmann–Jørgensen inequality to handle the squared summands. However, because the summands are not independent, we invoke Lemma 1. By applying this lemma, we obtain a partition \(\mathcal{G}\) of \(I_{{\boldsymbol{N}},{\boldsymbol{e}}}\) into \(|\mathcal{G}|=|I_{{\boldsymbol{N}},{\boldsymbol{e}}}|/n\) groups, each containing \(n\) i.i.d. observations. The i.i.d. property follows from the AHK representation 1 and the fact that within each group, any two observations share no common indices \(i_1,\dots,i_K\).

For each group \(G = \{{\boldsymbol{i}}_{1}(G), {\boldsymbol{i}}_{2}(G), \dots, {\boldsymbol{i}}_{n}(G)\} \in \mathcal{G}\), we define \[D_{f,{\boldsymbol{e}}}(G) = \frac{1}{n}\sum_{t=1}^n \Bigl(P_{{\boldsymbol{e}}}f\Bigr)^2\!\left(\{U_{{\boldsymbol{i}}_t(G)\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}\right),\] and let \[D_{f,{\boldsymbol{e}}} = \frac{1}{n}\sum_{t=1}^n \Bigl(P_{{\boldsymbol{e}}}f\Bigr)^2\!\left(\{U_{(t,\dots,t)\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}\right).\] Then we have \[\begin{align} \frac{1}{|I_{{\boldsymbol{N}},{\boldsymbol{e}}}|}\sum_{{\boldsymbol{i}}\in I_{{\boldsymbol{N}},{\boldsymbol{e}}}}(P_{{\boldsymbol{e}}}f )^2\left(\{U_{{\boldsymbol{i}}\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}\right)=\frac{1}{|I_{{\boldsymbol{N}},{\boldsymbol{e}}}|/n}\sum_{G\in \mathcal{G}} D_{f,{\boldsymbol{e}}}(G). \end{align}\] Note that for each \(G\in\mathcal{G}\), the AHK representation in 1 implies that \(D_{f,{\boldsymbol{e}}}\) and \(D_{f,{\boldsymbol{e}}}(G)\) are identically distributed. Consequently, by Jensen’s inequality, we have \[\begin{align} \mathbb{E}\Bigl[\sigma_{I_{{\boldsymbol{N}},{\boldsymbol{e}}}}^2\Bigr] &= \mathbb{E}\!\left[\left\|\frac{1}{|I_{{\boldsymbol{N}},{\boldsymbol{e}}}|/n}\sum_{G\in\mathcal{G}} D_{f,{\boldsymbol{e}}}(G)\right\|_\mathcal{F}\right] \le \mathbb{E}\!\left[\left\|D_{f,{\boldsymbol{e}}}\right\|_\mathcal{F}\right]. \end{align}\] Let us denote this bound by \[B_{n,{\boldsymbol{e}}} :=\mathbb{E}\!\left[\left\|D_{f,{\boldsymbol{e}}}\right\|_\mathcal{F}\right]= \mathbb{E}\!\left[\left\|\frac{1}{n}\sum_{t=1}^n \Bigl(P_{{\boldsymbol{e}}}f\Bigr)^2\!\Bigl(\{U_{(t,\dots,t)\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}\Bigr)\right\|_\mathcal{F}\right].\] Thus \(z\le\widetilde{z}:=\sqrt{B_{n,{\boldsymbol{e}}}}/\|(\pi_{\boldsymbol{e}})F\|_{P,2}\). Note that by symmetrisation inequality for independent processes, the contraction principle (Theorem 4.12. in [40]), and the Cauchy-Schwartz inequality, one has \[\begin{align} B_{n,{\boldsymbol{e}}}=&\mathbb{E}\left[\left\|\frac{1}{n}\sum_{t=1}^n (P_{{\boldsymbol{e}}}f )^2\left(\{U_{(t,...,t)\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}\right)\right\|_{\mathcal{F}}\right]\\ \le& \sigma_{\boldsymbol{e}}^2 +\mathbb{E}\left[\left\|\frac{1}{n}\sum_{t=1}^n\left\{ (P_{{\boldsymbol{e}}}f )^2\left(\{U_{(t,...,t)\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}\right)-\mathbb{E}\left[(P_{{\boldsymbol{e}}}f )^2\right]\right\}\right\|_{\mathcal{F}}\right]\\ \lesssim& \sigma_{\boldsymbol{e}}^2 +\mathbb{E}\left[\left\|\frac{1}{n}\sum_{t=1}^n\varepsilon_t\cdot (P_{{\boldsymbol{e}}}f )^2\left(\{U_{(t,...,t)\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}\right)\right\|_{\mathcal{F}}\right]\\ \lesssim& \sigma_{\boldsymbol{e}}^2 +\mathbb{E}\left[M_{\boldsymbol{e}}\left\|\frac{1}{n}\sum_{t=1}^n\varepsilon_t\cdot (P_{{\boldsymbol{e}}}f )\left(\{U_{(t,...,t)\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}\right)\right\|_{\mathcal{F}}\right]\\ \le& \sigma_{\boldsymbol{e}}^2 +\|M_{\boldsymbol{e}}\|_{P,2}\sqrt{\mathbb{E}\left[\left\|\frac{1}{n}\sum_{t=1}^n\varepsilon_t\cdot (P_{{\boldsymbol{e}}}f )\left(\{U_{(t,...,t)\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}\right)\right\|_{\mathcal{F}}^2\right]}. \end{align}\] An application of Hoffmann-Jørgensen’s inequality (Proposition A.1.6 in [36]) gives \[\begin{align} &\sqrt{\mathbb{E}\left[\left\|\frac{1}{n}\sum_{t=1}^n\varepsilon_t\cdot (P_{{\boldsymbol{e}}}f )\left(\{U_{(t,...,t)\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}\right)\right\|_{\mathcal{F}}^2\right]}\\ \lesssim& \mathbb{E}\left[\left\|\frac{1}{n}\sum_{t=1}^n\varepsilon_t\cdot (P_{{\boldsymbol{e}}}f )\left(\{U_{(t,...,t)\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}\right)\right\|_{\mathcal{F}}\right]+\frac{1}{n}\|M_{\boldsymbol{e}}\|_{P,2}. \end{align}\] By employing analogous reasoning to that used in the initial part of the proof, we deduce that \[\begin{align} &\mathbb{E}\left[\left\|\frac{1}{\sqrt{n}}\sum_{t=1}^n\varepsilon_t\cdot (P_{{\boldsymbol{e}}}f )\left(\{U_{(t,...,t)\odot {\boldsymbol{e}}'}\}_{{\boldsymbol{e}}'\le {\boldsymbol{e}}}\right)\right\|_{\mathcal{F}}\right]\\ \lesssim& \|P_{\boldsymbol{e}}F\|_{P,2}\int_{0}^{\widetilde{z}} \sup_{Q}\sqrt{1+\log N(P_{\boldsymbol{e}}\mathcal{F},\|\cdot\|_{Q,2},\epsilon \|P_{\boldsymbol{e}}F\|_{Q,2})}d\epsilon. \end{align}\] Note that the integral on the RHS can be bounded by \(J_{\boldsymbol{e}}(\widetilde{z})\) and thus \[\begin{align} B_{n,{\boldsymbol{e}}}\lesssim& \sigma_{\boldsymbol{e}}^2 + n^{-1}\|M_{\boldsymbol{e}}\|_{P,2}^2 + n^{-1/2} \|M_{\boldsymbol{e}}\|_{P,2} \|P_{\boldsymbol{e}}F\|_{P,2}J_{\boldsymbol{e}}(\widetilde{z}). \end{align}\]

Define \[\Delta=(\sigma_{\boldsymbol{e}}\vee n^{-1/2}\|M_{\boldsymbol{e}}\|_{P,2})/\|P_{\boldsymbol{e}}F\|_{P,2},\] it then follows that \[\begin{align} \widetilde{z}^2\lesssim \Delta^2 + \frac{\|M_{\boldsymbol{e}}\|_{P,2}}{\sqrt{n}\|P_{\boldsymbol{e}}F\|_{P,2}}J_{\boldsymbol{e}}(\widetilde{z}). \end{align}\] By applying Lemma 4 and Lemma 2.1 of [41] with \(J=J_{\boldsymbol{e}}\), \(A=\Delta\), \(B=\sqrt{\|M_{\boldsymbol{e}}\|_{P,2}/\sqrt{n}\|P_{\boldsymbol{e}}F\|_{P,2}}\) and \(r=1\), it yields that \[\begin{align} J_{\boldsymbol{e}}(z)\le J_{\boldsymbol{e}}(\widetilde{z})\lesssim J_{\boldsymbol{e}}(\Delta)\left\{1+J_{\boldsymbol{e}}(\Delta)\frac{\|M_{\boldsymbol{e}}\|_{P,2}}{\sqrt{n}\|P_{\boldsymbol{e}}F\|_{P,2} \Delta^2}\right\}. \end{align}\] Combining this with (16 ), we obtain the bound \[\begin{align} |I_{{\boldsymbol{N}},{\boldsymbol{e}}}|^{1/2}\mathbb{E}[\|H_{{\boldsymbol{N}}}^{\boldsymbol{e}}(f)\|_\mathcal{F}]\lesssim& J_{\boldsymbol{e}}(\Delta)\|P_{{\boldsymbol{e}}} F\|_{P,2} +\frac{J_{\boldsymbol{e}}^2(\Delta)\|M_{\boldsymbol{e}}\|_{P,2}}{\sqrt{n}\Delta^2}.\label{eq:local95maximal95ineq95main95bound} \end{align}\tag{17}\] Notice that \(\delta_{\boldsymbol{e}}\le \Delta\) by their definitions. By Lemma 4(iii), one has \[\begin{align} J_{\boldsymbol{e}}(\Delta)\le \Delta\frac{J_{\boldsymbol{e}}(\delta_{\boldsymbol{e}})}{\delta_{\boldsymbol{e}}}=\max\left\{J_{\boldsymbol{e}}(\delta_{\boldsymbol{e}}), \frac{\|M_{\boldsymbol{e}}\|_{P,2} J_{\boldsymbol{e}}(\delta_{\boldsymbol{e}})}{\sqrt{n}\|P_{\boldsymbol{e}}F\|_{P,2}\delta_{\boldsymbol{e}}}\right\}\le\max\left\{J_{\boldsymbol{e}}(\delta_{\boldsymbol{e}}), \frac{\|M_{\boldsymbol{e}}\|_{P,2} J_{\boldsymbol{e}}^2(\delta_{\boldsymbol{e}})}{\sqrt{n}\|P_{\boldsymbol{e}}F\|_{P,2}\delta_{\boldsymbol{e}}^2}\right\}, \end{align}\] where the second inequality follows from the fact \(J_{\boldsymbol{e}}(\delta_{\boldsymbol{e}})/\delta_{\boldsymbol{e}}\ge J_{\boldsymbol{e}}(1)\ge 1\). Finally, using Lemma 4(iii), \[\begin{align} \frac{J_{\boldsymbol{e}}^2(\Delta)\|M_{\boldsymbol{e}}\|_{P,2}}{\sqrt{n}\Delta^2}\le \frac{J_{\boldsymbol{e}}^2(\delta_k)\|M_{\boldsymbol{e}}\|_{P,2}}{\sqrt{n}\delta_k^2} \end{align}\] Combining the calculations with the bound in (17 ), we have the desired inequality. 0◻

Proof of Lemma 1↩︎

For \(K=1\), the result is trivial. For \(K\ge 2\), assume without loss of generality that \(N_1 \ge N_2 \ge \cdots \ge N_K.\) We prove for the case of \({\boldsymbol{e}}=(1,...,1)\) and \(I_{{\boldsymbol{N}},{\boldsymbol{e}}}=[{\boldsymbol{N}}]\) since other cases follow exactly the same arguments.

Our goal is to partition \([{\boldsymbol{N}}]\) into subsets (which we call groups) of size \(N_K\) that are transversal. For \(j=1,2,\dots,K-1\), define \(\phi_j\colon [N_K]\times [N_j] \to [N_j]\) by \[\phi_j(t,g)= \Bigl( ((t+g-2) \mod N_j) + 1\Bigr).\] That is, for each \(t\in [N_K]\) and \(g\in [N_j]\) the value \(\phi_j(t,g)\) is computed by adding \(t\) and \(g-1\), reducing modulo \(N_j\) (so that the result lies in \(\{0,1,\dots,N_j-1\}\)), and then adding 1 to get an element of \([N_j]\). For the \(K\)-th coordinate we set \(\phi_K(t) = t \quad \text{for } t\in [N_K].\) Index the groups by \[(g_1,g_2,\dots,g_{K-1})\in [N_1]\times [N_2]\times \cdots \times [N_{K-1}].\] Then, for each such \((g_1,\dots,g_{K-1})\), define \[G_{(g_1,\dots,g_{K-1})} = \Bigl\{\, \Bigl( \phi_1(t,g_1),\, \phi_2(t,g_2),\, \dots,\, \phi_{K-1}(t,g_{K-1}),\, t \Bigr) : t\in [N_K] \,\Bigr\}.\] Thus, each group contains \(N_K\) elements. We now claim that groups based on this mapping form a partition and satisfy transversality in each group.

We now claim the transversality property in each group. For a fixed group \(G_{(g_1,\dots,g_{K-1})}\) and a fixed coordinate \(j\) (with \(1\le j\le K-1\)), the \(j\)th coordinate of an element is given by \[\phi_j(t,g_j) = ((t+g_j-2) \mod N_j) + 1.\] Since the mapping \(t \mapsto ((t+g_j-2) \mod N_j) + 1\) is injective (note that \(N_K\le N_j\) so that there is no collision in the range), it follows that the \(j\)-th coordinates of the elements of \(G_{(g_1,\dots,g_{K-1})}\) are all distinct. For the \(K\)-th coordinate, the identity mapping \(t \mapsto t\) is trivially injective.

Next, we show the covering of \([{\boldsymbol{N}}]\) and disjointness of the groups. Recall that the total number of groups is \(N_1\cdot N_2\cdots N_{K-1}.\) Each group has \(N_K\) elements; hence, the union of all groups has \((N_1\cdot N_2\cdots N_{K-1})\cdot N_K = N\) elements. For surjectivity, let \(x = (x_1, x_2, \ldots, x_K)\) be an arbitrary element of \([{\boldsymbol{N}}] = [N_1] \times [N_2] \times \cdots \times [N_K]\). We wish to show that there exists \((g_1, g_2, \ldots, g_{K-1}) \in [N_1] \times [N_2] \times \cdots \times [N_{K-1}]\) and \(t \in [N_K]\) such that \[x = \Bigl( \phi_1(t, g_1),\, \phi_2(t, g_2),\, \dots,\, \phi_{K-1}(t, g_{K-1}),\, t \Bigr).\] Set \(t = x_K\). Then for each \(j = 1, 2, \ldots, K-1\), we must have \(\phi_j(x_K, g_j) = x_j,\) where by definition \(\phi_j(x_K, g_j) = (((x_K + g_j - 2) \bmod N_j) + 1)\). Notice that for each fixed \(x_K\), the mapping \[g \mapsto ((x_K + g - 2) \bmod N_j) + 1\] is an affine function (with coefficient \(1\)) on the cyclic group \(\mathbb{Z}/N_j\), and hence it is a bijection from \([N_j]\) onto \([N_j]\). Thus, for each \(j\) there exists a unique \(g_j \in [N_j]\) such that \(\phi_j(x_K, g_j) = x_j\). Therefore, every \(x \in S\) can be uniquely written in the form \[\Bigl( \phi_1(x_K, g_1),\, \phi_2(x_K, g_2),\, \dots,\, \phi_{K-1}(x_K, g_{K-1}),\, x_K \Bigr),\] which shows that the mapping \[\tau: (g_1,\dots,g_{K-1},t) \mapsto \Bigl( \phi_1(t, g_1),\, \phi_2(t, g_2),\, \dots,\, \phi_{K-1}(t, g_{K-1}),\, t \Bigr)\] is bijective.

Thus, the collection \[\mathcal{G} = \Bigl\{\, G_{(g_1,\dots,g_{K-1})} : (g_1,\dots,g_{K-1})\in [N_1]\times \cdots \times [N_{K-1}]\,\Bigr\}\] is a partition of \([{\boldsymbol{N}}]\) into groups of size \(N_K\), and in every group the entries in each coordinate are distinct. 0◻

6.6 Proof of Corollary 1↩︎

The proof follows the same arguments as in the proof of Corollary 5.3 in [14] using Lemma 5, and is not repeated here.

0◻

7 Auxiliary Lemmas↩︎

The following restates Theorem 5.2 in [22], which is a modification of Theorem 2.1 in [41] to allow for an unbounded envelope.

The following restates Lemma B.3 in [12].

The following is analogous to Lemma 5.2 in [14] and Lemma A.2 in [22].

The following restates a useful result of Lemma A.2. in [42] for the calculations of the uniform entropy integral of the classes after Hoeffding-type decomposition.

The following lemma is a multiway clustering version of Lemma 4.3 of [43].

The following lemma is a restatement of Lemma 3 in [6]. Here we consider the special case where each cell \(\boldsymbol{i} \in [\boldsymbol{N}]\) contains exactly one observation to avoid extra notations. Let \(I_1 = \bigcup_{\boldsymbol{e} \in \mathcal{E}_1} I_{\boldsymbol{N,e}}\) and \(k(\boldsymbol{i})\), denote the coordinate where \(\boldsymbol{i} \in I_1\) is non-zero.

References↩︎

---

1. First arXiv date: 11 Feb, 2026. We thank Jianfei Cao, Ulrich Hounyo, Michael Jansson, Michael Leung, and David Ritzwoller for their invaluable comments. This paper supersedes the earlier manuscript “Maximal inequalities for separately exchangeable empirical processes" (arXiv:2502.11432) by one of the authors. All remaining errors are our own.↩︎

2. Our framework accommodates high-dimensional regimes in which the array of data-generating processes may depend on the sample sizes; for notational economy, this dependence is left implicit.↩︎

3. See Section 2.3 in [16] for detailed explanation of empirical process notation for measure and integration.↩︎

4. For example, cluster-LASSO from [17] can be used for one-way clustering data. For two-way clustering panels, LASSO in [18] can be employed. For clustering more than two dimensions, the multiplier bootstrap for jointly exchangeable arrays in [12] can be used for choosing the valid penalty levels.↩︎

5. This condition need not hold uniformly in \(n\); in particular, \(\|F_n\|_{P,q}\to\infty\) is allowed.↩︎

6. See, e.g. [25].↩︎

7. Other candidates for inference procedures include the modified multiway empirical likelihood and the modified multiway jackknife variance estimator proposed in [30]. While these methods are computationally more demanding, they can potentially deliver improved higher-order asymptotic properties; see Theorem 3 therein.↩︎