Deep Neural Networks for Nonparametric Interaction Models with Diverging Dimension

DNN for Nonparametric Interaction Models

A]Sohom Bhattacharya[label=e1],
A]Jianqing Fan[label=e1]
A]Debarghya Mukherjee[label=e1]

[A]Department of Operations Research and Financial Engineering,
Princeton University[presep=, ]e1

Deep neural networks have achieved tremendous success due to their representation power and adaptation to low-dimensional structures.  Their potential for estimating structured regression functions has been recently established in the literature. However, most of the studies require the input dimension to be fixed and consequently ignore the effect of dimension on the rate of convergence and hamper their applications to modern big data with high dimensionality. In this paper, we bridge this gap by analyzing a $k^{th}$ order nonparametric interaction model in both growing dimension scenarios ($d$ grows with $n$ but at a slower rate) and in high dimension ($d \gtrsim n$).  In the latter case, sparsity assumptions and associated regularization are required in order to obtain optimal rates of convergence.
A new challenge in diverging dimension setting is in calculation mean-square error, the covariance terms among estimated additive components are an order of magnitude larger than those of the variances and they can deteriorate statistical properties without proper care.  We introduce a critical debiasing technique to amend the problem. We show that under certain standard assumptions, debiased deep neural networks achieve a minimax optimal rate both in terms of $(n, d)$. Our proof techniques rely crucially on a novel debiasing technique that makes the covariances of additive components negligible in the mean-square error calculation.  In addition, we establish the matching lower bounds.

Deep neural networks
High dimensional statistics
Non-parametric interaction model
Minimax rate
Sparse nonparametric components

§ INTRODUCTION

Recent advances in technology have allowed
statisticians to collect data on a large number of explanatory variables to
predict outcomes of interest [Goodfellow et al., 2016, Fan et al., 2020]. Often, the relationship between these outcomes and predictors are highly nonlinear (e.g., image data like MNIST [LeCun, 1998], CIFAR [Krizhevsky et al., 2009] etc.), yielding a practical need to investigate multivariate nonparametric regression model

\begin{equation}
Y_i= f(X_i)+ \varepsilon_i, \quad i=1,\ldots,n,
\end{equation}

where $X_i$'s are explanatory variables, $Y_i$'s are response variables and $\varepsilon_i$'s are unobserved errors. The statistical problem here is to recover $f$ under some minimal smoothness assumptions. A classical result of Charles Stone [Stone, 1982] shows if $f$ is a $d$-variate $(\beta,C)$ smooth function (precise definition given later), the minimax optimal rate of estimation is of the order $n^{-\frac{2\beta}{2\beta+d}}$, which is referred to as “curse of dimensionality", i.e., large sample size is necessary to estimate the regression function well. In particular, when $\beta$ is fixed, it is easily to see that when $d \asymp \log n$, minimax consistent estimators can not be obtained.

To tackle aforementioned issue, one needs to impose a low-dimensional structure on the function $f$ [Kpotufe and Dasgupta, 2012, Yang and Dunson, 2016, Yang and Tokdar, 2015], such as additive regression model [Stone, 1985], projection pursuit regression model [Friedman and Stuetzle, 1981], higher order interaction model [Stone, 1994], generalized higher-order interaction model [Horowitz and Mammen, 2007, Kohler and Krzyżak, 2016], single and multi-index models [Hardle et al., 1993], etc.; see [Fan and Gijbels, 1996] for an overview. The key idea behind all these simplifications is to reduce the complexity of the underlying function class and mollify the effect of dimension in the rate of convergence. For example, [Stone, 1985] proved that for additive regression model (when $f(x) = \sum_{j=1}^d f_j(x_j)$) with $\beta$-smooth univariate components, the minimax optimal rate of estimation of $f$ in terms of squared $L_2(P_X)$ loss is $C_dn^{-\frac{2\beta}{2\beta + 1}}$. [Fan et al., 1998] demonstrated further the adaptivity of such an additive structure and [Horowitz and Mammen, 2007] studied further a general class of nonparametric regression models with unknown link functions. As each component function $f_j$ is univariate and can be estimated at a rate $n^{-\frac{2\beta}{2\beta + 1}}$, the effect of dimension appears only through the multiplicative constant $C_d$, not in the power of $n$. Later, [Stone, 1994] extended this result for higher-order interaction model (i.e. when $f(x) = \sum_{J \in S} f_J(x_J)$, $S$ is a collection of subsets of $\{1, \dots, d\}$ and $|J| \le d^\star$ for all $J \in S$) and showed that  the minimax optimal rate of estimation is $C_dn^{-\frac{2\beta}{2\beta + d^\star}}$, where, as before, the effect of $d$ appears through a multiplicative constant. Again, all of the above results assume that $d$ is finite, which is inappropriate in many modern data science applications.

In practice, it is not enough just to obtain minimax optimal rates, efficient and easily computable estimators are warranted. Several methods like kernel, spline, and wavelet based estimators (see Section 2 of [Fan and Gijbels, 1996] for an overview), backfitting algorithm [Buja et al., 1989],  off-the-shelf nonparametric methods like boosting [Freund et al., 1996], random forest [Breiman, 2001] have been popularly used. Recently, neural networks have emerged as an imperative tool for analyzing nonlinear relations between covariates and response variables. Deep (or multilayer) neural networks have been the backbone of the incredible advances in machine learning that have resulted in massive success in reinforcement
learning, robotics, computer
vision, natural language processing, and other statistical prediction tasks <cit.>.
Leveraging the availability of large amounts of digitized data, deep learning has enjoyed a plethora of empirical successes. However, most of the successes are purely human engineered (i.e. achieved after tuning lots of hyperparameters), lacking theoretical guarantees. As a consequence, researchers are naturally interested in understanding the statistical properties of such useful and practical estimators. In [Kohler and Krzyżak, 2005], the authors established that for the standard nonparametric models as well as for the higher-order interaction models, deep neural network (henceforth DNN) based estimators are minimax rate optimal up to logarithmic factor. Later, the results are extended for a more general class of functions (namely generalized hierarchical interaction model, see [Bauer and Kohler, 2019]) in a series of work by[Kohler and Krzyżak, 2016, Bauer and Kohler, 2019, Schmidt-Hieber, 2020, Kohler and Langer, 2021, Fan et al., 2022].

The previous research in understanding the rate of convergence of DNN-based estimators assumes the underlying dimension of covariate $d$ is fixed and consequently ignores the effect of $d$ in the rate of the estimator.

However, in the era of big data, the dimension of underlying covariates is quite large in many practical problems, often larger than the underlying sample size $n$. For example, in a simple image classification problem, a $28 \times 28$ image lies in $\reals^{784}$. Similarly, in genome-wide association studies, the number of single nucleotide polymorphism (SNP) can be significantly higher than the number of individuals[Novembre et al., 2008].  In that paper, the authors have analyzed a dataset consisting of $1387$ individuals and $\sim 2\times 10^5$ SNPs. In such examples,

a sharper analysis quantifying the optimal dependence on $d$ is often necessary.

In this paper, we aim to bridge this gap by analyzing the nonparametric regression model both in growing dimension (when $d$ grows with $n$ but $d = o(n)$) and high dimension (when $d \gtrsim n$) setup. The distinction between these two cases lies in whether or not the sparsity and regularization are required for consistent estimation, in an analogue way to the linear model with growing dimension and sparse linear model in high dimension.  We analyze the following $k$-way interaction model [Stone, 1994]:

\begin{equation}\label{eq:define_f}
Y_i = \sum_{l=1}^k \sum_{\substack{(j_1 < \dots <j_l) \in [d]^l}}f_{j_1,j_2,\ldots,j_l}(X_{i,\{j_1,j_2,\ldots,j_l\}}) + \eps_i \,,
\end{equation}

where $[d]^l$ is the collection of all ordered subsequence of length $l$ of $\{1, \dots, d\}$.  As mentioned before, $k = o(\log n)$ is necessary in order to have a consistent estimator.  Hence, we will take a finite $k$.
Such models include the well-studied additive model [Stone, 1985], namely, $f(x)=\sum_{j=1}^d f_j(x_j)$ and interaction models [Stone, 1994]. When $d \gtrsim n$, it is necessary to impose sparsity assumption even in the high-dimensional linear model. We, therefore, consider a sparse $k$-way interaction model as follows:

\begin{equation}\label{eq:define_spam}
Y_i = \sum_{l=1}^k \sum_{\substack{(j_1,j_2,\ldots,j_l) \in S_l}}f_{j_1,j_2,\ldots,j_l}(X_{i,\{j_1,j_2,\ldots,j_l\}}) + \eps_i \,,
\end{equation}

where $S_l \subset [d]^l$ with $|S_l| \ll d^l$.

When each univariate function is constrained to be linear, the model (<ref>) reduces to a sparse parametric interaction model. When $k=1$, the model (<ref>) becomes a sparse additive model which has been well-studied in literature [Lin and Zhang, 2006, Koltchinskii and Yuan, 2010, Ravikumar et al., 2009, Tan and Zhang, 2019, Yuan and Zhou, 2016], building upon the recent developments on penalized linear regression. For example, [Tan and Zhang, 2019] proposes to doubly penalize each component $f_j$ by its empirical norm and functional semi-norm to estimate the regression function.

When $d \ll n$, we estimate the mean function of (<ref>) via minimizing squared error loss (details can be found in Section <ref>). However, when $d \gtrsim n$, we need to penalize explicitly to enforce sparsity. In previous works on the sparse additive models, researchers typically impose two penalties: i) one to enforce sparsity (e.g. via the summation of $L_2(\bbP_n)$ norm on the component functions) and ii) another to control the complexity of the component functions (e.g., penalize with respect to RKHS norm if the component functions lie in a RKHS). In reality, implementing such a doubly penalized estimator is often computationally challenging. To overcome such difficulty, we introduce a two-step hard thresholding-based estimator, which is motivated by the seminal work of SURE independent screening [Fan and Lv, 2008, Fan and Song, 2010] and least square estimation after model selection in high dimensional linear regression model [Belloni and Chernozhukov, 2013]. The idea can be briefly described as follows: first, we obtain an initial estimator (say $\hat f^{\text{init}}$) of the mean function by only penalizing the summation of empirical $L_2(\bbP_n)$ norm of the component functions (which can be implemented by penalizing the last layer of neural network). Analogously, this can also be thought of as a version of a group lasso penalty, where the neural network corresponding to each component is a group of variables and we penalize the sum of the norm of each group (here we take the $L_2(\bbP_n))$ norm) to enforce a component-wise sparsity.
Next, we further perform a hard thresholding on the non-zero component of $\hat f^{\text{init}}$ to weed out the small non-zero components and estimate the active sets. Finally, we perform an empirical risk minimization by minimizing the squared error loss over only the components selected in the previous step to obtain the final estimator. We prove that under some fairly standard assumptions, the estimator is minimax rate optimal. We now summarize our contribution as follows.

Contribution:  Our main contribution is a rigorous theoretical analysis of $k$-way interaction model (also known as nonparametric ANOVA model) in both growing and high dimensional setup. To achieve this goal, we also prove several results that may be of independent interest. We summarise the key contributions below:

* We analyze the nonparametric $k$-way interaction model using neural network based estimator both when $d = \dim(X)$ increases with $n$ and $d = o(n)$, and when $d \gtrsim n$ with aforementioned regularization. We show that the neural network-based estimator is minimax rate optimal up to a poly-log factor.

* We introduce a novel debiasing technique that makes the covariances among additive components negligible, which reduces statistical errors.

* We prove an approximation result of smooth function using deep neural network (Theorem <ref>) via a novel debiasing technique, which implies that one can approximate a smooth function using neural networks at the same rate even under constraints (i.e. marginals of a multivariate function are $0$).

* We establish the minimax lower bound for estimating the high dimensional $k$-way interaction model, which, to the best of our knowledge, is not present in the literature.

The rest of the paper is organized as follows: In Section <ref>, we analyze the $k$ way interaction model when $d \ll n $. We divide the entire analysis into three subsections: Subsection <ref> bounds approximation error, Subsection <ref> analyzes the statistical error, and Subsection <ref> develops the minimax lower bound. Section <ref> deals with the analysis when $d \gtrsim n$. We broadly divide our analysis into two parts; Subsection <ref> contains the analysis for fixed designs and Subsection <ref> deals with the case of random designs. Furthermore, in Subsection <ref>, we establish the minimax lower bound for the sparse nonparametric interaction models.Section <ref> provides conclusional remarks. Finally, Section <ref> contains the proof of Theorem <ref>. The remaining proofs can be found in the Appendix.

§ ANALYSIS OF DNN IN LOW DIVERGING DIMENSION

In this section, we present our analysis of the DNN-based estimator of the mean function when the dimension $d$ of $X$ grows with $n$ but $d \ll n$. We consider a $k$-order interaction nonparametric regression model. For an input-output pair $(X, Y)$, where $X \in \reals^d$ and $Y \in \reals$, a $k$-order interaction model is defined as:

\begin{equation}
\label{eq:model_k}
Y_i = \sum_{\substack{J \subset [d] \\ |J| \le k}}f_{0, J}(X_{i,J}) + \eps_i \,.
\end{equation}

Here, $\eps_i$'s are assumed to be centered error independent of $X$ and the functions $f_{0, J}: \reals^{|J|} \mapsto \reals$. For simplicity of exposition, we henceforth confine ourselves to a $2$-order interaction model. The extension of our analysis from a $2$-order interaction model to a general $k$-order interaction model is purely technical and the proof ideas will be outlined at the end of the relevant sections. For $k = 2$, the the model presented in equation (<ref>) simplifies to:

\begin{equation}
\label{eq:two_way_ss_anova}
Y_i = \sum_{j=1}^d f_{0, j}(X_{i,j}) + \sum_{j < k}f_{0, jk}(X_{i, j}, X_{i, k}) + \eps_i \equiv f_0(X_i) + \varepsilon_i \,.
\end{equation}

where $f_{0, j}'s$ are univariate functions and $f_{0, jk}$'s are bivariate functions. We will make the following assumption on $X$ and $\varepsilon$:

$X$ is supported on $[0, 1]^d$ and admits a density function $p$ such that $\sup_{x \in [0,1]^d}p(x)=: p_{\max} < \infty$, where $p_{\max}$ is free of $d$. We assume $\eps$ is sub-gaussian with sub-gaussian constant $\sigma^2_\eps$.

Assumption <ref> is a standard assumption in the literature of nonparametric regression. Note that assuming $X$ has compact support is equivalent to assuming that the support is $[0, 1]^d$ via centering and scaling. Many of our results can be extended to the case of unbounded $X$ via truncation arguments – we omit such arguments for the sake of the simplicity of our article. The upper bound assumption on the density is quite natural, as this merely rules out certain degenerate distributions (i.e. when the density of some region diverges to infinity).  The assumption of sub-gaussianity on the error terms is natural as we aim to analyze the least square estimate of $f$ based on neural networks. It is well-established in the literature that for heavy-tailed errors, a least square estimator of a nonparametric function is not minimax optimal (e.g., see [Han and Wellner, 2019]) and one should use a variant of Huber loss function. We believe most of our analysis can be extended to this heavy-tailed error, but this is out of the scope of the current paper. The recent article [Fan et al., 2022] sheds some light on the estimation error of deep neural networks with heavy-tailed error in fixed dimension setup.

In our model (<ref>), the component functions of $f_0$ are not identifiable without further assumptions for two reasons:

* All the functions are identifiable up to shift, i.e. we cannot identify the difference between $(f_{0, i_1 i_2}, f_{0, j_1 j_2})$ and $(f_{0, i_1 i_2} +c, f_{0, j_1 j_2} - c)$ for a constant $c$.

* The univariate marginals of the bivariate functions are not identifiable. For example, consider the subcollection of functions $\cC_1 = \{f_{0, 1}, f_{0, 12}, f_{0, 13}, \dots, f_{0, 1d}\}$. Let $g_k(x) = \int f_{0, 1k}(x, y) \ dy$ be the marginal of $f_{0, 1k}$. Then we cannot differentiate between $\cC_1$ and a modified collection $\cC_2 = \{f_{0, 1} + \sum_{k \geq 2}g_k, f_{0, 12} - g_2, \dots, f_{0, 1d} - g_d\}$.

Therefore, to identify and estimate the non-parametric functions, we need to impose certain structural conditions which will prevent us from shifting. Toward that end, we assume the following:
[Identifiability and boundedness]

We assume the following conditions for identifiability purpose and boundedness:

* All the univariate functions in (<ref>) integrates to $0$, i.e. $\int_0^1 f_{0, j}(x) \ dx = 0$ for $1 \le j \le d$.

* All the bivariate functions have $0$ marginals, i.e.

\int_0^1 f_{0, ij}(x, y) \ dx = \int^1_0 f_{0, ij}(x, y) \ dy = 0 \,.

* $\|f_0\|_\infty \le B$ for some $B > 0$, where $f_0$ is defined in (<ref>).

Note that, the above assumption implies that the total integral of $f_{0, ij}$ is also $0$. Assumption <ref> ensures identifiability of all the functions involved in (<ref>). Our next assumption is a smoothness assumption on the underlying component function:

The functions $\{f_{0, j}\}$ and $\{f_{0, ij}\}$ are assumed to be $(\beta, L)$-Holder smooth, i.e. the functions are $\lfloor \beta
\rfloor$ times differentiable and $\lfloor \beta \rfloor^{th}$ derivative are Holder with exponent $\beta - \lfloor \beta \rfloor$ and constant $L$.

The smoothness assumption is standard in the nonparametric regression literature, cf. <cit.>, as this smoothness assumption controls the complexity of the underlying function class. Another common assumption is that all $f_{0,i}$ and $f_{0,ij}$ belong to a Reproducing Kernel Hilbert Space(RKHS), e.g, see [Raskutti et al., 2012], [Koltchinskii and Yuan, 2010] and references therein. It would be interesting to study the analog of our results under such assumptions.

Assuming that $\{f_{0, ij}\}$ have a different level of smoothness $\beta_{ij}$ is equivalent to assuming they all have $\min_{i,j} \beta_{i,j}$ smoothness as long as one is concerned about estimation error bounds.

To estimate the mean function of (<ref>) using deep neural networks, we need some properties of DNNs. A neural network is a parametric function $f_\theta$ which maps the input space to output space, i.e. in our context $f_\theta$ maps $\reals^d$ to $\reals$. A two-layer (i.e. one hidden layer) neural network (often termed as shallow neural network) with $N$ neurons and activation function $\sigma$ is defined as:

f_\theta(x) = a + \sum_{i=1}^N b_i \sigma(c_i^\top x + d_i) \,.

where $\theta = \{a, b_i, c_i, d_i\}$. Therefore, $f_\theta$ first projects $x$ into $\reals^N$ via a linear transformation $x \mapsto Cx + d$ where $c_i^\top$'s are the rows of $C$ and $d_i$'s are elements of $d$. Then it applies a non-linear activation $\sigma$ to all the coordinates of $Cx + d$. Finally, it takes a linear combination of the coordinates using the map $x \mapsto a + \langle b, x
\rangle$. This shallow neural network can be extended to the deep neural network by adding more hidden layers. A graphical representation of the flow of neural networks is depicted in Figure <ref>.

shorten >=1pt,->,
node distance=,
every pin edge/.style=<-,shorten <=1pt,
neuron/.style=circle,fill=black!25,minimum size=17pt,inner sep=0pt,
input neuron/.style=neuron, fill=blue!50,
output neuron/.style=neuron, fill=blue!50,
hidden neuron/.style=neuron, fill=green!50,
annot/.style=text width=4em, text centered

/ in 1,...,1

[input neuron, pin=left:Input] (I-) at (0,-1.5) ;

in 1,...,
in 1,...,3
node[hidden neuron] (H-) at (*,-cm) ;

[annot,above of=H-1, node distance=1cm] (hl) Hidden layer ;

/ in 1,...,1

[output neuron, pin=right:Output] (O-) at (4*,-1.5) ;

in 1,...,1
in 1,...,3
(I-) edge (H1-);

[remember=as (initially 1)] in 2,...,
in 1,...,3
in 1,...,3
(H-) edge (H-);

[remember=as (initially 3)] in ,...,
in 1,...,3
in 1,...,1
(H-) edge (O-);

[annot,left of=hl1] Input layer;
[annot,right of=hl] Output layer;

A graphic representation of a deep neural network with 3 hidden layers, one input layer, one output layer.

For a neural network, we denote by $L$ the number of hidden layers (termed as depth) and by $N$, the maximum number of neurons at the hidden layers (denoted by width). Henceforth, we denote by $\cF_{NN}(d, N, L, W, o)$ by the collection of all neural networks with depth $L$, width $N$, the total number of active (non-zero) weights $W$, input dimension $d$ and output dimension $o$. We may sometimes omit the input and output dimensions in the specification of the class of neural networks when it is clear from the context.

The expressibility of neural networks has been a topic of interest for decades, which started to flourish at the beginning of 1990s. It was proved in [Hornik et al., 1989] that the set of all shallow neural networks (i.e. with one hidden layer) is dense in the space of all Borel measurable functions and any Borel measurable function $f$ can be estimated within arbitrary accuracy by increasing the number of neurons in the hidden layer. The degree of accuracy was quantified by [Barron, 1993] – if the Fourier transform of the function has a finite first moment, then the approximation error decreases at a rate $1/\sqrt N$, where $N$ is the number of hidden neurons. The error of approximating any smooth function by a multi-layer feed-forward network was established by [Mhaskar, 1996] along with a choice of $N, L$ which achieves the optimal accuracy. Since then, there has been a series of research on the approximation capability of deep neural networks in terms of their width, height, weights, and activation function – recently, researchers have obtained precisely bound on the approximation error of smooth functions in terms of depth and width of neural network  cf. [Yarotsky, 2017, Lu et al., 2021, Fan and Gu, 2022] and references therein. For example, Theorem 1.1 of [Lu et al., 2021] shows that if $f:[0, 1]^d \to \reals$ is $\beta$ times differentiable with bounded derivatives then there exists a neural network $\phi$ with width $c_1 N\log{N}$ and depth $c_2 L\log{L} + d$ ($c_1, c_2$ are some constants depending on $\beta$) with the following approximation error:

\|f - \phi\|_\infty \le C(NL)^{-\frac{2\beta}{d}} \,,

where the constant $C$ depends on $(d, \beta)$.

A related line of research delves into understanding how DNNs can successfully adapt to the unknown underlying low dimensional structure of the mean function in nonparametric regression. Performing nonparametric regression using suitable DNNs and achieving minimax optimal error bound (up to some polylog factor) initially started from the pioneering work of [Kohler and Krzyżak, 2005]. The minimax rate of estimation of a mean function $f$ in a standard nonparametric regression problem (e.g., additive noise model $Y = f(X) + \eps$) is $n^{-2\beta/(2\beta +d)}$ where $\beta$ is the smoothness index of $f$ (Assumption <ref>) and $d$ is the dimension of $X$. The estimation error suffers from the curse of dimensionality, i.e. the rate is very slow for large $d$. However, as mentioned in the introduction, it is possible to circumvent this curse by imposing more structure on $f$ such as  additive (see [Stone, 1985]), or some higher order interaction model (see [Stone, 1994]).

In [Kohler and Krzyżak, 2005], the authors show that it is possible to adapt to the rate $n^{-2\beta/(2\beta + k)}$ using neural network-based estimate if the model is a $k$-way interaction model. The key advantage of using a neural network is that we only need to specify its architecture, i.e. the width and layer, not the exact structure of the mean function and the estimated NN becomes minimax rate optimal. Recently, in a series of papers [Kohler and Krzyżak, 2016, Schmidt-Hieber, 2020, Bauer and Kohler, 2019, Kohler and Langer, 2021, Fan et al., 2022, Fan and Gu, 2022], the authors also establish that it is possible to estimate the underlying mean function $f_0$ at a rate $n^{-2\beta/(2\beta + k)}$ via neural networks when $f_0$ belongs to a more general class than $k$-way interaction model, called generalized hierarchical interaction model.

All the above analyses assume that the underlying dimension of $d$ is fixed. This is at odds with many modern applications.  Often time the constant in front of the rate of estimation depends on the dimension of $X$ even for the simple additive model as we need to estimate $d$ univariate functions. When $d$ is fixed, it is possible to get away with a cruder bound. However, if $d \rightarrow \infty$, we need new techniques to obtain the optimal dependence on $d$, which is  one of the primary goals of this paper.

A schematic diagram for the structured deep neural networks that are used to estimate the structured nonparametric regression with interactions. The first part consists of additive fully connected neural networks for the univariate predictors and the second component comprises the summation of fully connected bivariate neural networks for approximating bivariate interaction effects.

We now briefly describe our estimation procedure. Consider the model in equation (<ref>). We select two neural networks $(\hat \phi_1, \hat \phi_2)$ via minimizing the squared-error loss:

\begin{equation}\label{eq:estimator_growing_d}
\left(\hat \phi_1, \hat \phi_2\right) = \argmin_{\phi_1 \in \cF_{NN}^1,\phi_2 \in \cF_{NN}^2}\frac1n \sum_{i=1}^n \left(Y_i - \phi_1(X_i) - \phi_2(X_i)\right)^2
\end{equation}

and set estimate $\hat f_{\text{big}} := \hat \phi_1 + \hat \phi_2$. Finally set the estimator $\hat f : = sgn (\hat f_{\text{big}})(|\hat f_{\text{big}}| \wedge B)$, i.e. truncate the output of $\hat f$ at $[-B, B]$ as we assume $\|f_0\|_\infty \le B$ (see point 3 of Assumption <ref>). Here the class of neural networks $\cF^1_{NN}$ and $\cF^2_{NN}$ are defined as:

\begin{align}
\label{eq:def_F1} \cF_{NN}^1 & = \cF_{NN}\left(d, c_1 dN_1\log{N_1}, c_2L_1\log{L_1}, c_3 d(N_1\log{N_1})^2L_1\log{L_1}, 1\right) \\
\label{eq:def_F2} \cF_{NN}^2 & = \cF_{NN}\left(d, c_1 d^2N_2\log{N_2}, c_2L_2\log{L_2}, c_3 d^2(N_2\log{N_2})^2L_2\log{L_2}, 1\right)
\end{align}

for some constants $c_1, c_2, c_3$. For example, in $\cF^2_{NN}$, all pairwise interactions are used as input variables.  Note that we are not using a fully connected neural network here. Rather, we estimate each univariate (resp. bivariate) component function via fully-connected neural networks of width $N_1 \log{N_1}$ (resp. $N_2\log{N_2}$) and depth $L_1 \log{L_1}$ (resp. $L_2\log{L_2}$) and then add them.
See Figure <ref> for an illustration.  With slight abuse of notation, we still use $\cF_{NN}^1$ and  $\cF_{NN}^2$ to denote these specially structured neural networks.
The quantities $N_i$ and $L_i$ typically depend on the sample size $n$ and will be specified later. Like any other non-parametric estimator, the generalization error of $\hat f$ can be decomposed into two parts:

\begin{equation}\label{eq:error_decomp}
\|\hat f - f_0\|_{L_2(P_X)} \le \underbrace{\|\hat f - f^\star\|_{L_2(P_X)}}_{\text{Statistical error}}  + \underbrace{\|f^\star - f_0\|_{L_2(P_X)}}_{\text{Approximation error}}
\end{equation}

Here $f^\star = \phi_1^\star + \phi_2^\star$ is (approximately) the best approximator of $f_0$ in terms of $L_2(P_X)$ norm among the class of neural networks, i.e.

\|f^\star - f_0\|_{L_2(P_X)} \le \inf_{\phi_1 \in \cF_{NN}^1, \phi_2 \in \cF_{NN}^2}\|f_0 - \phi_1 - \phi_2\|_{L_2(P_X)} + \tau_n

for some small tolerance $\tau_n$ to be specified later.
The rest of our analysis is divided into two parts: Subsection <ref> bounds the approximation error and Subsection <ref> bounds the statistical error in terms of $(d, N_i, L_i, n)$. Finally, in Subsection <ref> we choose the architectures of $\cF_{NN}^1, \cF_{NN}^2$ to balance the errors. This will prove that the rate obtained by the neural network is minimax optimal both in terms of $(n, d)$, up to some logarithmic factor.

It is important to note, we are estimating $f_0$ as the sum of two neural networks with different architectures instead of using a single neural network. To see why this is necessary, note that the rate of convergence of the estimator is $dn^{-2\beta_1/(2\beta_1 + 1)}$ and $d^2n^{-\beta_2/(\beta_2 + 1)}$ (Corollary <ref>) upto a logarithmic factor, where $\beta_1$ is the smoothness index of the univariate component and $\beta_2$ is the smoothness index of the bivariate components of $f_0$. Now note that if $\beta_2$ is sufficiently large, then $dn^{-2\beta_1/(2\beta_1 + 1)}$ will be dominant rate, otherwise $d^2n^{-\beta_2/(\beta_2 + 1)}$ will be the dominant one. If we only use one neural network, then the rate obtained is $d^2(n^{-2\beta_1/(2\beta_1 + 1)}  + n^{-\beta_2/(\beta_2 + 1)})$, which is not minimax optimal (Theorem <ref>) especially when $\beta_2 \gg \beta_1$. The trade-off between the error terms is more complicated as $d:=d(n) \rightarrow \infty$ and using two neural networks (one for univariate components and the other for bivariate components) we can obtain the optimal error bounds.
We conjecture that it is not possible to obtain the minimax error bound by using a single neural network.

§.§ Approximation error

In this section, we compute the approximation error of mean estimation in $2$-way interaction model (<ref>).

To this end, we invoke Theorem 1.1 of [Lu et al., 2021], which states any $\beta$ smooth function $f:[0, 1]^d \to \reals$ can approximate $f$ within error of $(NL)^{-2\beta/d}$ in $L_\infty$ norm using a DNN with width $O(N\log{N})$ and depth $O(L\log{L} + d)$. Recall that, by Assumption, <ref>, $\int f_j(x)dx=0$ and both the marginal integrals of $f_{ij}$s are $0$. As stated in the previous section, we will approximate the $f_i$s and $f_{ij}$s by two separate neural networks with different depths and widths. Before stating our main theorem, we first quantify a bound on the growth of $d$ in comparison to $n$ in the following assumption:

We assume that $d$ is growing with $n$ and

\left(dn^{-\frac{2\beta_1}{2\beta_1 + 1}} + \binom{d}{2}n^{-\frac{\beta_2}{\beta_2 + 1}}\right)\log^{4.5}{n} = o(1) \,.

The following theorem establishes the approximation error of the mean function via neural networks:

Consider the two-way interaction model in the equation
where all the univariate components are $\beta_1$ smooth and all the bivariate components are $\beta_2$ smooth.  Choose $N_1 L_1 = \lfloor n^{1/2(2\beta_1 + 1)}\rfloor$, $N_2 L_2 = \lfloor n^{1/2(\beta_2 + 1)}\rfloor$ Then we have:

\begin{align*}
& \inf_{\phi_2 \in \cF^1_{NN}, \phi_2 \in \cF^2_{NN}}\bbE\left[\|f_0(X) - \phi_1(X) - \phi_2(X)\|^2\right] \\
& \qquad \qquad \qquad \le C_3 \left(d(N_1L_1)^{-4\beta_1} + \binom{d}{2} (N_2L_2)^{-2\beta_2}\right) \,,
\end{align*}

where $\cF^1_{NN}$ and $\cF^2_{NN}$ are same as defined in (<ref>) and (<ref>).
All the constants $(C_1, C_2, C_3)$ are independent of $(d, N_1, L_1, N_2, L_2)$.

The proof of Theorem <ref> is deferred to Section <ref>. Here we provide the sketch of the proof. The key idea in our proof is to construct a neural network that not only approximates $f$ but also (approximately) satisfies the identifiability constraint Assumption <ref>, i.e. the marginals of $f_{ij}$'s are $0$. This will nullify the effect of the higher-order terms yielding the optimal rate of approximation. For univariate components, this integral constraint can be satisfied exactly; consider a function $f:[0, 1] \to \reals$ where $f$ satisfies Assumption <ref> with some $\beta$ and $\int_0^1 f(x) \ dx = 0$. Then by <cit.>, there exists a neural network $\phi$ with width of $O(N_1\log{N_1})$ and depth $O(L_1\log{L_1})$ such that
$$\|f - \phi\|_\infty \le C(N_1L_1)^{-2\beta}$$
where the constant $C$ depends only on $\beta$ and Sobolev norm of $f$.
Define $I_\phi := \int_0^1 \phi(x) \ dx$ and define a new neural network $\tilde \phi = \phi - I_\phi$ guaranteeing $\int \tilde \phi=0$. As subtracting a constant from a neural network amounts to changing the bias of the last layer, the architecture of $\phi$ and $\tilde \phi$ are the same. Using triangle inequality,

\|f - \tilde \phi\|_\infty \le 2C(N_1L_1)^{-2\beta} \,.$$

The above idea doesn't work for the approximation of the bivariate functions $f_{ij}$ as we not only need the integral to be $0$ but also the marginals to be $0$, i.e. $\int f_{ij}(x, y) \ dx = \int f_{ij}(x, y) \ dy = 0$. To address this issue, we devise a debiasing technique summarized as follows: given any bivariate function $f$ which is $\beta$ times differentiable with bounded derivatives, there exists a neural network $\phi(x, y)$ (using <cit.>) with architecture $O(N_2\log{N_2}), O(L_2\log{L_2})$ that satisfies:

\|f - \phi\|_\infty \le C(N_2L_2)^{-\beta} \,.

Given such a neural network $\phi$, we construct $\phi_1$ (resp. $\phi_2$) such that $\phi_1(x)$ (resp. $\phi_2(y)$) is a neural network with width $O(N_2\log{N_2})$ and depth $O(L_2\log{L_2})$
which approximates $\int \phi(x, y) \ dy$ (resp. $\int \phi(x, y) \ dx$). Our final estimator is of the form $\tilde \phi = \phi - \phi_1 - \phi_2 + \iint \phi(x, y) \ dx \ dy$. We show that $\|\tilde \phi- f\|_{\infty}=O((NL)^{-\beta})$ and $\int \tilde \phi(x) \ dx = \int \tilde \phi(x, y) \ dy = O\left((NL)^{-2\beta}\right)$.
This smaller approximation error ($(N_2L_2)^{-2\beta}$ instead of $(N_2L_2)^{-\beta}$) is crucial to obtain the optimal error bounds. Details of the proof can be found in Appendix <ref>.

we want to investigate the effect of dimension on the rate of convergence of the deep neural network in a nonparametric regression problem,  $Y_i = f_0(X_i) + \eps_{i}$, $i \in [n]$ where $Y_i \in \mathbf{R}, X_i \in \mathbf{R}^d$ and the errors $\varepsilon_i$ are i.i.d. sub-Gaussian with parameter $\sigma^2$. We work in the regime $d=o(n)$. We will need to further assume an upper bound on the density of $X_i$, given as follows.

The goal is to determine the optimal rate of convergence with respect to the neural network classifier.

Towards that goal, we want to construct a class of deep neural network estimators $\mathcal{F}_n$ such that $d(\hat{f}_n,f_0)$ is small, where $d$ is the $L_2(P_X)$ distance and $f_n$, is defined by

\begin{equation}\label{eq:fn_define}
\hat f_n = \argmin_{f \in \cF_n} \frac1n \sum_{i=1}^{n} (Y_i - f(X_i))^2 \,.
\end{equation}

Given the set $\cF_n$, we will provide an upper bound on $d(\hat{f}_n,f_0)$ using the following Theorem whose proof we include in the Appendix for the helpfulness of the readers.

Consider a  non-parametric linear regression model $Y = f_0(X) + \eps$. Suppose we estimate $f_0$ using the Sieve type of technique, i.e. at stage $n$, we find the least squared estimator over a class of functions $\cF_n$ given by (<ref>).
Further define the approximation error is $\omega_n$, i.e. $\inf_{f \in \cF_n} d(f, f_0) \le \omega_n$ and $V_n$ is the VC-dimension of $\cF_n$. Assume that all the functions in $\cF_n$ and $f_0$ are uniformly bounded by $1$. Then we have:

d(\hat f_n, f_0) = O_p\left(\sqrt{\frac{V_n}{n}\log{\frac{n}{V_n} }\sqrt{\log{n}}} + 2\omega_n\right)

upto log-factor, where $d$ is the $L_2(P_X)$ distance.

§ NOTATIONS

For a function $f:\mathbb{R}^d \mapsto \mathbb{R}$, define $\|f\|^2_n= \frac{1}{n} \sum_{i=1}^{n}f^2(X_i)$, and $\|f\|^2_2= \int f^2 dP_X$. Given two sequences $a_n,b_n \ge 0$, we say $a_n \lesssim b_n$ if $\exists C>0$ such that $a_n \le C b_n$.
Let $\alpha=q+s$ for some $q \in \mathbb{N}_0$ and $0 \le s<1$. A function $f:\mathbb{R} \mapsto \bbR$ will be called $(\alpha,C)$-smooth, if for every $\boldsymbol\alpha=(\alpha_1,\ldots,\alpha_d)\in \mathbb{N}^d_0$ with $\sum_{j=1}^{d}=q$, the partial derivatives $tial^q f/ tial x^{\alpha_1}_1 \ldots tial x^{\alpha_d}_d$ exist and

\[
\left| \frac{tial^q f}{tial x^{\alpha_1}_1 \ldots tial x^{\alpha_d}_d} (x)- \frac{tial^q f}{tial x^{\alpha_1}_1 \ldots tial x^{\alpha_d}_d} (z)\right| \le C \|x-z\|^s
\]

§.§ Statistical error and main theorem

In this section, we discuss the techniques to bound the statistical error of estimation and also present the main theorem that bounds the out-of-sample prediction error. Statistical error of an estimator typically relies on the number of samples and the complexity of the underlying class of function. Recall from (<ref>) that we approximate $f_0$ by the sum of two neural networks $\hat \phi_1 + \hat \phi_2$ where $\hat \phi_1$ approximates the univariate component and $\hat \phi_2$ approximates the bivariate component. It is immediate from the proof of Theorem <ref> that, the approximation error depends on $N$ and $L$ through their product $NL$. Therefore, we work with the neural network with constant depth and optimize over width. So, $\cF^1_{NN}$ and $\cF^2_{NN}$ (equation (<ref>) and (<ref>)) can be revised as:

\begin{align}
\label{eq:modified_def_F1} \cF_{NN}^1 & = \cF_{NN}\left(d, c_1 dN_1\log{N_1}, c_2, c_3 d(N_1\log{N_1})^2, 1\right), \\
\label{eq:modified_def_F2} \cF_{NN}^2 & = \cF_{NN}\left(d, c_1 \binom{d}{2}N_2\log{N_2}, c_2, c_3 \binom{d}{2}(N_2\log{N_2})^2, 1\right).
\end{align}

In practice, one estimates $\hat \phi_1, \hat \phi_2$ by optimizing over the weights of the neural network through gradient descent (e.g. see 6.2 of [Goodfellow et al., 2016]). Since our problem is highly non-convex in terms of the weights of the neural network, it is not immediately clear whether such an algorithm will converge to global minima. To avoid such issues, we work with a neural network that minimizes empirical risk in our article. Understanding the properties of the gradient descent-based estimator is left open for future research.

The aim of this section is to bound the variance term, i.e. $\bbE[(\hat f(X) - \phi^\star(X))^2 \mid \bS_n]$ where $\phi^\star$ is the best approximator of the mean function $f_0$ and $\bS_n:= \{(X_i,Y_i), i \in [n]\}$ be the training sample. The prediction/test error depends on the trade-off between the complexity of the underlying function class and the number of training samples. If the underlying function class is too complex compared to the number of samples $n$, then it may overfit the data resulting in high generalization error. Therefore, it is imperative to quantify the complexity of the class of neural networks. The following lemma of [Bartlett et al., 2019] establishes a bound on the VC dimension of the class of neural networks with the number of active weights $W$ and depth $L$:

Suppose $V(W,L)$ denotes the largest VC- dimension of a ReLU network with $W$ parameters and $L$ layers. There exist constants $c_1, c_2>0$ such that

\[ c_1 WL \log(WL) \le V(W,L) \le c_2 WL \log{WL}. \]

Roughly speaking, if we have a fully connected neural network with width $N$ and depth $L$, then $W \sim N^2L$ and consequently the VC dimension is $\sim N^2L^2$.
Once we have a bound on the complexity of the underlying function class, we can use techniques from empirical process theory to obtain the rate of convergence of the generalization error. For a VC-type function class with VC dimension $V$, the statistical error is $O(\sqrt{V/n})$ up to logarithmic factors.
To see the necessity of the assumption, consider the rate obtained in Corollary <ref>. Assumption <ref> ensures the consistency of the estimator. On the other hand, it is also necessary up to the logarithmic factor as we establish in Section 4 that the rate obtained in Corollary <ref> is minimax optimal up to the logarithmic factor (Theorem <ref>). Therefore the assumption is the weakest possible assumption on the growth of $d$ with respect to $n$ under no further structural assumption. The following theorem presents a bound on the overall approximation error, combining both statistical and approximation errors.

Under Assumption <ref> - <ref> and <ref>, the ERM estimator $\hat f$ of $f_0$ defined in equation (<ref>) satisfies:

\bbE\left[\left(\hat f(X) - f_0(X)\right)^2 \mid \bS_n\right] = O_p\left(\rho_n^2 + \frac{V_n}{n}\log^{3/2}{n}\right),
%O_p\left(\left(dn^{-\frac{2\beta_1}{2\beta_1 + 1}} + \binom{d}{2}n^{-\frac{\beta_2}{\beta_2 + 1}}\right)\log^{4.5}{n}\right) \,.

\begin{align*}
\rho_n^2 &= \text{approximation error} \le C_1\left(dN_1^{-4\beta_1} + \binom{d}{2}N_2^{-2\beta_2}\right) \\
V_n &= \text{complexity} \le C_2\left(dN_1^2\log^2{N_1}\log{(dN_1)} + \binom{d}{2}N_2^2\log^2{N_2}\log{(dN_2)}\right) \,.
\end{align*}

The proof of the theorem can be found in Appendix <ref>. Now we need to choose the value $N_1, N_2$ carefully to balance both the approximation error and statistical error. This leads to the following corollary:

Choosing $N_1 = n^{1/2(2\beta_1+ 1)}$ and $N_2 =n^{1/2(\beta_2 + 1)}$ (take the nearest integer if they are not integers) in Theorem <ref>, we have:

\bbE\left[\left(\hat f(X) - f_0(X)\right)^2 \mid \bS_n\right] =O_p\left(\left(dn^{-\frac{2\beta_1}{2\beta_1 + 1}} + \binom{d}{2}n^{-\frac{\beta_2}{\beta_2 + 1}}\right)\log^{4.5}{n}\right) \,.

In the subsequent section, we show that this rate is minimax optimal up to the logarithmic factor. An immediate question that stems from this analysis is whether we can improve the dependence on the logarithmic factor. The recent article [Fan and Gu, 2022] proposes a slightly different approximation technique than that of [Lu et al., 2021], which might reduce the power of the logarithmic factor from $4.5$.
However, the more interesting question is whether we can completely get rid of the log factor which is still open.

[Extension for k-way interaction model]
Our analysis for the two-way interaction model can also be extended verbatim to the $k$-way interaction model. Here we provide a sketch of the extension. For identifiability of the model, we need to assume all marginals of $j$-component functions are $0$, $j \ge 2$. We estimate $f_0$ via the sum of $k$ neural networks $\phi_1 + \dots + \phi_k$ where

\phi_i \in \cF_{NN}\left(d, c_1\binom{d}{i}N_i\log{N_i}, c_2, \binom{d}{i}(N_i\log{N_i})^2, 1\right) \,.

Since the debiasing technique used in our proofs rely on the approximation of polynomials by neural networks, they can be extended to a general $k$-way interaction model. This in turn controls the approximation error. The statistical error can be bounded similarly by the VC-dimension of the neural networks via Lemma <ref>.

§.§ Minimax lower bound

In this section, we establish the minimax lower bound for estimating the non-parametric $2$-way interaction model (<ref>), when the dimension of the covariate $d= d(n) \rightarrow \infty, d=o(n)$ and satisfies Assumption <ref>.

For our estimation problem, the minimax risk is defined as:

\mathfrak{M}(n, d, \cF) = \inf_{\hat f} \sup_{\substack{f \in \cF \\ X \sim P_{X}}}\bbE_f\left[\left(\hat f(X) - f(X)\right)^2\right]

Here $f(X) = \sum_{j=1}^d f_j(X_j) + \sum_{i < j} f_{ij}(X_i, X_j)$ and the supremum is taken over all $\{f_j\}$ and $\{f_{ij}\}$ where the component functions belongs $\Sigma(\beta, L)$ (see Assumption <ref>).

The main difficulty in establishing the minimax lower bound is to incorporate the effect of growing dimension $d$. Our problem is certainly harder than the problem of estimating $f_1(x_1) + f_{12}(x_1, x_2)$, assuming all other components are known in advance.  The issue is how other components contribute to the lower bound in the growing $d$ setting.  The following theorem gives a precise answer.

Consider the two-way interaction model as defined in equation (<ref>) where each component functions $\{f_i\}_{i=1}^d$ and $\{f_{ij}\}_{1 \le i < j \le d}$ belongs to $\Sigma(\beta, L)$ (see Assumption <ref>) and denote this collection of mean functions as $\cF$. Then the minimax rate of estimation under Assumptions <ref>-<ref> is:

\mathfrak{M}(n, d, \cF) = \inf_{\hat f} \sup_{\substack{f \in \cF \\ X \sim P_{X}}}\bbE_f\left[\left(\hat f(X) - f(X)\right)^2\right] \ge c\left( dn^{-\frac{2\beta_1}{2\beta_1 + 1}} +  \dbinom{d}{2}n^{-\frac{2\beta_2}{2\beta_2 + 2}}\right) \ \,.

for some constant $C$ independent of $(n, d)$.

The proof of this theorem can be found in Appendix <ref>. A few remarks are in order; first of all, Theorem <ref> establishes the fact that neural network-based estimate of the mean function $f_0$ is minimax rate optimal up to a poly-log factor. This result complements the one in [Kohler and Langer, 2021, Schmidt-Hieber, 2020] for the non-parametric regression in a fixed dimension regime. Secondly, Theorem <ref> is derived under Assumption <ref>, where we assume all univariate components are $\beta_1$-smooth and all bivariate components are $\beta_2$-smooth.
However, our proof can be easily adapted to the setting where each function has different smoothness and we have to pay the price for the lowest smoothness. See Remark <ref> where we have pointed out the key steps to prove such extension.

[Extension for k-way interaction model]
The proof of the minimax lower bound for the two-way interaction model can be easily generalized for the general $k-$way interaction model by constructing alternatives similarly and using Fano's inequality. The optimal rate is given by

\inf_{\hat f} \sup_{f \in \cF, X \sim P_{X}} \bbE_{P_{X}}[(\hat f(X) - f(X))^2] \gtrsim \left(\sum_{j=0}^k \dbinom{d}{j}n^{-\frac{2\beta_j}{2\beta_j + j}}\right)

We have studied the effect of dimensionality in the minimax upper and lower bound for the interaction model in our article. A natural next step is to find Pinsker's constant for the $k$-way interaction model in diverging dimensions, which is left as a future research direction.

§ ANALYSIS OF DNN WHEN $D \GG N$

This section presents our analysis when the underlying dimension of $X$ is larger than the sample size. Although simpler parametric models have been well-studied to avoid the curse of dimensionality, the literature on nonparametric estimation for $k$-way interaction model (<ref>) in this regime is relatively sparse. In the nonparametric framework, a significant amount of research has been done on a high dimensional additive model under sparsity constraint where the model under consideration is as follows:

Y_i = \sum_{j \in S}f_{0, j}(X_{ij}) + \eps_i

where $S \subseteq \{1, \dots, d\}$ with $s = |S| \ll n$. This is a standard assumption in high dimensional statistical analysis, where the true signal depends on a few covariates, but the active set $S$ is apriori unknown. The analysis of the sparse additive model in high dimension was initiated in [Lin and Zhang, 2006] and later refined in a series of papers literature [Koltchinskii and Yuan, 2010, Ravikumar et al., 2009, Tan and Zhang, 2019, Yuan and Zhou, 2016]. [Koltchinskii and Yuan, 2010] obtained error bounds under a global boundedness assumption, which was subsequently removed by [Raskutti et al., 2012]. Recently, [Tan and Zhang, 2019] obtained minimax guarantees when the underlying component functions lie in a reproducing kernel Hilbert space. All the previous works typically use two penalties for optimal estimation: one to control the complexity of the underlying function class and the other to control the sparsity. Yet, it is enough to use only one penalty for inducing sparsity for a neural network-based estimator, as the complexity can be controlled through the network's width and depth.

In spite of such extensive works on sparse linear models, the study of the high-dimensional linear models with $k$-way interaction ($k>1$) is very sparse, letting alone $k$-way non-parametric interactions.   For some references about the variable selection in high dimensional linear interaction models, readers may consult [Zhao et al., 2009], [Bien et al., 2013], [Hao and Zhang, 2014] and references therein. In this section, we bridge this gap by studying the asymptotic properties of the neural network estimator.
As before, we only elaborate on the analysis of the $2$-way interaction model and comment on how to extend our analysis for the general $k$-way interaction model. The two-way sparse additive model is defined as follows:

\begin{equation}\label{eq:define_sparse_model}
Y_i= \sum_{j \in S_1}f_{0, j}(X_{ij})+ \sum_{(k,l)\in S_2} f_{0, kl}(X_{ik},X_{il}) + \varepsilon_i \equiv f_0(X_i)+\varepsilon_i,
\end{equation}

where $S_1 \subset [d]$, $S_2 \subset [d]\times[d]$ with $S_1,S_2$ sparse, i.e., $s_i:=|S_i| \ll n$ for $i=1,2$. Define this class of functions by $\mathcal{F}_{\text{sp}}$. Here also we assume that the univariate components are $\beta_1$ smooth and the bivariate components are $\beta_2$ smooth. From Theorem <ref>, there exists $\{\phi^\star_j\}_{j \in S_1}, \left\{\phi^\star_{kl}\right\}_{(j, k) \in S_2}$, $j\in S_1, (j,k) \in S_2$, such that, $\phi^\star_j \in \mathcal{F}_{NN,1}$, $\phi^\star_{jk} \in \mathcal{F}_{NN,2}$ and

\begin{align*}
& \|f_{0,j} - \phi^\star_{j}\|_\infty \le C_1 N_1^{-2\beta_1} \ \forall \ j \in S_1 \,,  \\
& \|f_{0,kl} - \phi^\star_{kl}\|_\infty \le C_2 N_2^{-\beta_2} \ \forall \ (k < l) \in S_2 \,.
\end{align*}

\begin{align}
\label{eq:def_nn_uni}\cF_{NN, 1} & = \cF_{NN}(1, c_1N_1\log{N_1}, c_2, c_3 N^2_1\log^2{N_1}, 1) \\
\label{eq:def_nn_bi} \cF_{NN, 2} & = \cF_{NN}(2, c_4 N_2\log{N_2}, c_5, c_6 N^2_2\log^2{N_2}, 1) \,.
\end{align}

Finally, define

\begin{align}\label{eq:nn_approx_sparse}
\phi^\star=\sum_{j \in S_1}\phi^\star_j +\sum_{(k<l) \in S_2} \phi^\star_{kl}.
\end{align}

Following Assumption <ref>, we assume the following growth condition on the sparsity parameter $s_1, s_2$:

The sparsity parameters $s_1, s_2$ corresponding to the univariate components and bivariate components satisfy the following condition:

s_1\left(n^{-\frac{2\beta_1}{2\beta_1 + 1}}\log^4{n} + \frac{\log{d}}{n}\right) + s_2\left(n^{-\frac{\beta_2}{\beta_2 + 1}}\log^4{n} + \frac{\log{d}}{n}\right) = o(1) \,.

The reason behind Assumption <ref> is similar for that of Assumption <ref>, i.e. without this assumption, there will be no consistent estimator (modulo the logarithmic factor) as it is minimax optimal rate of estimation. The analysis for the high dimensional model is not a straightforward extension of the techniques used when $d = o(n)$. To understand why, recall that even in a simple sparse linear model with $\ell_1$ penalty, it is not possible to obtain a minimax rate optimal estimator without a form of restricted strong convexity assumption, which is typically not needed when the dimension grows slower than the sample size. Therefore, for ease of presentation, we divide the entire analysis into three parts: Section <ref> establishes the rate of convergence for the fixed design model, i.e. when $X_i$'s are some fixed points in $\reals^p$. Section <ref> deals with the random design, i.e. $X$ is assumed to be a random variable. Finally, in Section <ref> we present the minimax lower bound to establish that our neural network-based estimator is minimax rate optimal up to log factors.

§.§ Analysis of Fixed design model

In this section, we present our analysis for the fixed design model, where we assume $X_i$'s to be fixed. Similar to (<ref>), we estimate univariate and bivariate components using neural networks with different architecture, but at the same time, we use the $\ell_1$-norm of $L_2(\bbP_n)$ penalty to enforce sparsity.
See page 966 of [Antoniadis and Fan, 2001] and [Yuan and Lin, 2006] for such an idea in selecting a group of variables.
Our estimator $\hat f = \sum_j \hat \phi_j + \sum_{k < l} \hat \phi_{kl}$ where components are defined as

\begin{align}
\label{eq:fixed_sparse_penalize}
\left\{\hat \phi_j\right\}, \left\{\hat \phi_{kl}\right\} & = \mathop{\arg \min}\limits_{\substack{\phi_j \in \cF_{NN, 1} \\ \phi_{kl} \in \cF_{NN, 2}}} \left[\frac1n \sum_i \left(Y_i - \sum_j \phi_j(X_{ij}) - \sum_{k < l}\phi_{kl}(X_{ik}, X_{il})\right)^2 \right.\nonumber \\
& \qquad \qquad \qquad \qquad \qquad \left. + \left(\lambda_{n,1}\sum_j \|\phi_j\|_n + \lambda_{n,2}\sum_{k < l} \|\phi_{kl}\|_n\right)\right]
\end{align}

where $\cF_{NN, 1}$ and $\cF_{NN, 2}$ are defined in (<ref>) and (<ref>) respectively, $\| \cdot \|_n$ is the $L_2(\bbP_n)$ norm and $\lambda_{n,1}, \lambda_{n,2}$ to be specified later.We now make couple of remarks on the estimation procedure. First, although we put $L_2(\bbP_n)$ penalty on the component functions to enforce sparsity, we may use the $\ell_1$ penalty on the last layer of the weights (which adds the component neural networks). Such a penalty, in spite of its computational benefit, is difficult to analyze theoretically. Second, as we do not necessarily assume $\beta_1=\beta_2$, we need two regularization parameters $\lambda_{n, 1}$ and $ \lambda_{n, 2}$. We also need different architectures to adapt to the smoothness as highlighted by Remark <ref>. Furthermore, for technical simplicity, we assume the following $l^\infty$ bound:

\begin{equation}\label{eq:B_bound}
\|f_j\|_{\infty} \le B, \qquad \|f_{jk}\|_{\infty}\le B.
\end{equation}

Consequently, for numerical stability, we truncate the component neural networks (i.e. $\phi_j, \phi_{jk}$) at level $B$.
The following theorem establishes the rate of convergence.

Under Assumption <ref>-<ref> and <ref> along with (<ref>),
the estimator obtained in (<ref>)  satisfes

\|\hat f - f_0\|_n^2 = O_p\left(s_1\left(\rho_{n, 1}^2 +\lambda_{n, 1} + \frac{\lambda_{n,1}^2}{2}\right)  + s_2 \left(\rho_{n, 2}^2  + \lambda_{n,2} + \frac{\lambda_{n, 2}^2}{2}\right)\right)

where $\rho_{n, 1}$ (resp. $\rho_{n, 2}$) is the approximation error of the univariate (resp. bivariate) components of the mean function by neural networks, bounded by (<ref>), provided that

\begin{equation}\label{eq:define_lambda}
\lambda_{n, 1} = C_3\sqrt{\frac{V_{n, 1} \log{n}}{n} + \frac{2\log{d}}{n}}, \quad
\lambda_{n, 2} = C_4\sqrt{\frac{V_{n, 2} \log{n}}{n} + \frac{3\log{d}}{n}}   \,,
\end{equation}

with $V_{n, 1}  = N_1^2 \log^3{N_1}$ and $V_{n, 2} = N_2^2 \log^3{N_2}$.

Let us try to understand and simplify the rate obtained by Theorem <ref>. Recall from Theorem <ref>, we have:

\begin{equation} \label{eq:approx-err}
\rho_{n, 1} \le C_1 N_1^{-2\beta_1}, \quad \rho_{n, 2} \le C_2 N_2^{-\beta_2} \,,
\end{equation}

for some constants $C_1, C_2>0$.

Furthermore, it is revealed in our proof (see  (<ref>) and (<ref>)) that  $\lambda_{n, 1}$ and $\lambda_{n, 2}$ given by (<ref>) are the optimal choices, in which $V_{n, 1}$ and
$V_{n, 2}$ are the order of VC-dimensions for $\cF_{NN, 1}$ and $\cF_{NN, 2}$ respectively.
Now, typically $\lambda_{n, i} = o(1)$ (which holds as soon as $V_{n,i} = o(n/\log n)$, a condition necessary for consistency) and consequently $\lambda_{n,i}^2=o(\lambda_{n, i})$. Hence the rate of convergence in Theorem <ref> can be simplified as

\|\hat f - f_0\|_n^2 = O_p\left(s_1\left(\rho_{n, 1}^2 +\lambda_{n, 1} \right)  + s_2 \left(\rho_{n, 2}^2  + \lambda_{n,2}\right)\right) \,.

Now we choose $N_1, N_2$ to balance $\rho_{n, i}^2$ and $\lambda_{n, i}$, which leads to the following corollary:

Choosing $N_1 = n^{1/(2(1 + 4\beta_1))}$ and $N_2 = n^{1/(2(1 + 2\beta_2))}$ (take the nearest integer if they are not integers), we obtain:

\|\hat f - f_0\|_n^2 = O_p\left(s_1\left(n^{-\frac{2\beta_1}{1 + 4\beta_1}}\log^2{n} + \sqrt{\frac{\log{d}}{n}}\right)  + s_2 \left(n^{-\frac{\beta_2}{1 + 2\beta_2}}\log^2{n} + \sqrt{\frac{\log{d}}{n}}\right)\right) \,.

cAs we will see in Subsection <ref>, this rate is not minimax optimal. The reason is similar to that for the high dimensional sparse regression model: if we do not assume any condition on the curvature of the loss function (e.g. restricted isometry or restricted eigenvalue type conditions) then it is not possible to obtain any minimax optimal estimator which is computable in polynomial time as proved in [Zhang et al., 2014].
Hence, we need to assume a form of Restricted Strong Convexity (RSC) to obtain a minimax optimal error bound. RSC-type assumptions were popularized by [Bickel et al., 2010, Candès and Tao, 2007]. A similar assumption is also used in the analysis of the sparse additive regression model, e.g. see <cit.>. In particular, we assume the following:

There exist constants $\kappa_1, \kappa_2 > 0$, such that for any function $\phi = \sum_j \phi_j +\sum_{k<l} \phi_{kl}$
that satisfies:

\begin{align}
& \lambda_{n, 1} \sum_{j \in S_1^c}\|\phi_j - \phi^\star_j \|_n + \lambda_{n, 2}\sum_{(k < l) \in S^c_2}\|\phi_{kl} - \phi^\star_{kl}\|_{n} \notag \\
& \le 4 (s_1\rho_{n, 1}^2 + s_2 \rho_{n, 2}^2)  + 3\lambda_{n, 1} \sum_{j \in S_1}\|\phi_j - \phi^\star_j \|_{n} +  3\lambda_{n, 2} \sum_{(k < l) \in S_2}\|\phi_{kl} - \phi^\star_{kl}\|_{n} \notag \\
\label{eq:cond_for_rsc} & \hspace{5em} + s_1\lambda_{n, 1}^2 + s_2\lambda_{n, 2}^2
\end{align}

also satisfies:

\begin{equation}
\label{eq:main_rsc_condition}
\kappa_1^2\sum_{j \in S_1}\|\phi_j - \phi^\star_j\|^2_{n} + \kappa_2^2 \sum_{(k < l) \in S_2}\|\phi_{kl} - \phi^\star_{kl}\|^2_{n} \le \|\phi - \phi^\star\|_n^2
\end{equation}

where $\phi^\star$ is defined by (<ref>), $\rho_{n, 1}$ and $\rho_{n,2}$ are defined in Theorem <ref> and $\lambda_{n, 1}$ and $\lambda_{n,2}$ are defined in (<ref>).

Before going into further details, we first compare our assumption with the standard restricted eigenvalue (RE) condition for the high dimensional linear model [Bickel et al., 2009]. Roughly speaking, the RE condition assumes:

\begin{equation}
\label{eq:re_explanation}
\sum_{j \in S} (\beta_j - \beta_{0, j})^2 \lesssim \frac1n \|X(\beta - \beta_0)\|_2^2
\end{equation}

for any $\beta$ that satisfies:

\begin{equation}
\label{eq:cone_cond}
\sum_{j \in S^c} |\beta_j| \lesssim \sum_{j \in S}|\beta_j - \beta_{0, j}| \,,
\end{equation}

where $S$ is the true active set. Now consider the nonparametric additive model, where $X_j$ affects $Y$ through $f_j(X_j)$. Furthermore, we typically use the sum of $L_2(\bbP_n)$ norm as a penalty function in place of $L_1$ norm for sparsity. Ignoring the approximation error for the time being (i.e. assuming true $f_j$ belong to the function class over which we optimize), a natural analog of condition (<ref>) is the following:

\begin{equation}
\label{eq:re_explanation_nonparam}
\sum_{j \in S} \|f_j - f_{0, j}\|_n^2
\lesssim \|\hat f - f_0\|_n^2
\end{equation}

for all $f = \sum_j f_j$ that satisfies (similar to (<ref>))

\begin{equation}
\label{eq:cone_nonparam}
\sum_{j \in S^c}\|f_j \|_n \lesssim \sum_{j \in S} \|f_j - f_{0, j}\|_n \,.
\end{equation}

This is precisely what was used in <cit.> for analyzing additive models in high dimensions. Assumption <ref> is almost the same as this one with some modification due to the fact we are approximating the component functions through neural networks and consequently we may have non-zero approximation error. Ignoring the bivariate components, (<ref>) is exactly the same as (<ref>) with $f$ replaced by $\phi$ and $f_0$ replaced by $\phi^\star$, which is the best approximator of $f_0$ from the class of neural networks. The only difference is in
(<ref>), which differs from (<ref>) as we need to take into account the approximation error. This is precisely why the term $s_1 \rho_{n, 1}^2 + s_2 \rho_{n, 2}^2 $ appears on the right-hand side of equation (<ref>). The other additional term $s_1 \lambda_{n, 1}^2 + s_2 \lambda_{n, 2}^2$ is due to some mathematical artifact of the proof and also can be made arbitrarily smaller in order. For simplicity, if we ignore the bivariate components and this additional term for the time being, then (<ref>) simplifies to:

\begin{align}
\label{eq:cond_for_rsc_additive}
\lambda_{n, 1} \sum_{j \in S_1^c}\|\phi_j - \phi^\star_j \|_n  & \le 4 s_1\rho_{n, 1}^2 + 3\lambda_{n, 1} \sum_{j \in S_1}\|\phi_j - \phi^\star_j \|_{n} \,.
\end{align}

If the function class is well-specified (as in the case of linear regression or additive models in [Tan and Zhang, 2019]), then the first term of the RHS can be removed. Consequently, cancelling $\lambda_{n, 1}$ from both sides, our condition becomes same as (<ref>) as used in [Tan and Zhang, 2019].

We now state the statistical error theorem:

Consider the function $\hat{\phi}$ defined by (<ref>). Then, under the same assumptions as that of Theorem <ref> along with Assumption <ref>, we have:

\begin{equation}
% \frac{1}{n}\sum_{i=1}^{n} (f_0(X_i) -\hat{\phi}(X_i))^2 \lesssim s_1 ((NL)^{-4\beta}+\lambda^2_{n,1})+ s_2 ((NL)^{-2\beta}+\lambda^2_{n,2}).
\|\hat f - f_0\|_n^2 = O_p\left(s_1 (\rho_{n, 1}^2+\lambda^2_{n,1})+ s_2 (\rho_{n, 2}^2+\lambda^2_{n,2})\right) \,.
\end{equation}

The proof of this theorem is presented in  Appendix <ref>. The key difference between the rate obtained in Theorem <ref> and Theorem <ref> is the absence of $\lambda_{n, i}$, which yields a faster rate under Assumption <ref>. Similar discussion as that of after Theorem <ref> is in order; we need to choose $N_1, N_2$ to balance $\rho_{n, i}^2$ and $\lambda_{n, i}^2$. This leads to the following corollary:

Choosing $N_1 = n^{1/(2(1 + 2\beta_1))}$ and $N_2 = n^{1/(2(1 + \beta_2))}$ (take the nearest integer if they are not integers), we obtain:

\|\hat f - f_0\|_n^2 = O_p\left(s_1\left(n^{-\frac{2\beta_1}{1 + 2\beta_1}}\log^4{n} + \frac{\log{d}}{n}\right)  + s_2 \left(n^{-\frac{\beta_2}{1 + \beta_2}}\log^4{n} + \frac{\log{d}}{n}\right)\right) \,.

In Section <ref>, we prove a matching lower bound to obtain that this rate is in fact minimax optimal up to log factors.

§.§ Random design setting

In the previous section, we have analyzed the fixed design model, i.e. $X_i$'s are assumed to be fixed. Here, we extend our analysis to the random design model, namely when $X_i$'s are random variables with distribution satisfying Assumption <ref>.

As $X_i$'s are assumed to be random, we analyze the behavior of an estimator in terms of expected squared error loss, also known as generalization/out-of-sample error. Let us first highlight the difference between this subsection and the previous one. In the fixed design setup, the entire analysis hinges on the given $(X_1, \dots, X_n)$, i.e. we don't require to evaluate the performance of the predictor on some unobserved $X$. However, when we assume $X$'s are random, then it is imperative to evaluate the performance on the unseen $X$'s as there is a high probability (in fact probability is 1 if some component of $X$ is continuous) that in future we need to predict on new observations. Hence, in this setup we need to bound the expected squared error loss, i.e. $\|\hat f - f_0\|_{L_2(P_X)}^2$ instead of empirical squared error loss, i.e. $\|\hat f - f_0\|_n^2$. A key step to bound this generalization error is to control the fluctuation of the corresponding empirical process, which ensures that a predictor, that performs well on the training samples, also performs well on the previously unobserved test samples. Typically, it is difficult to achieve such a guarantee using only the $L_1$-norm of the $L_2(\bbP_n)$ penalties as this penalty only controls the complexity of the functions on the observed samples. Although it might be possible to obtain a minimax optimal estimator by using $L_\infty$ penalty along with the $L_1$-norm of the $L_2(\bbP_n)$ penalties on the component functions, the optimization procedure becomes computationally challenging. To circumvent this issue, we propose a computationally efficient two-step procedure, which, at the cost of mildly stronger assumptions, allows us to estimate the true mean function optimally. We illustrate our estimator via Algorithm <ref>.
Constant $c_1, c_2$ and penalty parameters $\lambda_{n, 1}$, $\lambda_{n, 2}$.
$\hat f$: the estimator based on neural network.
Dataset $\{(X_1, Y_1), \dots, (X_n, Y_n)\}$.

Divide the dataset into two equal halves with $n/2$ data in each set (if $n$ is odd, then take $(n+1)/2$ data in the first half and $(n-1)/2$ in the second half). Denote the halves by $\cD_1$ and $\cD_2$.

Using $\cD_1$, estimate the component functions $\{\hat \phi^{\text{init}}_j\}$ and $\{\hat \phi^{\text{init}}_{jk}\}$ by solving (<ref>).

Set $\hat S_1 = \{j: \|\hat \phi^{\text{init}}_j\|_n \ge c_1\lambda_{n, 1}\}$ and $\hat S_2 = \{(j < k): \|\hat \phi^{\text{init}}_{jk}\|_n \ge c_2 \lambda_{n, 2}\}$.

Re-estimate component functions on the estimated active set $\hat S_1$ and $\hat S_2$ by minimizing (un-penalized) squared error loss:

\begin{align}
\label{eq:two_step_OLS_eq}
& \{\hat \phi^{\text{final}}_j\}_{j \in \hat S_1}, \{\hat \phi^{\text{final}}_{jk}\}_{(j, k) \in \hat S_2} \notag \\
& = \argmin_{\substack{\phi_j \in \cF_{NN, 1} \\ \phi_{j, k} \in \cF_{NN,2}}} \frac{1}{n}\sum_i \left(Y_i - \sum_{j \in \hat S_1}\phi_j(X_{ij}) - \sum _{(j, k) \in \hat S_2}\phi_{jk}(X_{ij}, X_{ik})\right)^2 \,.
\end{align}

Return the final estimate $\hat f^{\text{final}} = \sum_{j \in \hat S_1}\hat \phi^{\text{final}}_j + \sum_{(j < k) \in \hat S_2}\hat \phi^{\text{final}}_{jk}$  .
Estimation under random design setting

The key idea is as follows: we first split the samples into two (almost) equal halves. Based on the first half of the sample, we estimate the component functions by the same penalized procedure used in the fixed design setting (see(<ref>)). Next, we use hard thresholding on the estimators of univariate and bivariate component functions to estimate the active set (Step 3 of Algorithm <ref>). Denote the estimated active sets by $\hat S_1$ and $\hat S_2$ respectively. Note that our threshold levels are proportional to the penalty applied in (<ref>). The constants $c_1, c_2$ will be mentioned explicitly in the proof. Once we have $\hat S_1, \hat S_2$, then we solve the unpenalized least square problem only on the active set based on the second half of the data (Step 4 of Algorithm <ref>) to obtain the final estimate. Our algorithm is primarily motivated by SURE independent screening  [Fan and Lv, 2008, Fan and Song, 2010] and the idea of least square regression upon selecting an active set via LASSO or other penalized regression  [Belloni and Chernozhukov, 2013]. SURE independent screening selects the active variables by performing marginal regression for
additive models, whereas for sparse linear regression, [Belloni and Chernozhukov, 2013] proposes to choose the active set using LASSO and then perform OLS on the selected subset to reduce bias. We combine these two ideas in Algorithm <ref>. First, we estimate the active subset of univariate and bivariate components $\hat S_1, \hat S_2$ (step 2-3 of Algorithm <ref>) and ensure that with high probability: i) $S_i \subseteq \hat S_i$ and ii) $|\hat S_i \cap S_i^c| = O(|S_i|)$. We need i) to avoid any potential bias that may occur by ignoring active components and ii) to guard against false positives, i.e. we do not select too many inactive variables. To ensure this, we need the following assumption:

Assume that all the component functions has minimal signal strength $r_n$, i.e. $\min_{j \in \cS_1} \|f_j^0\|_2 \ge r_n$, $\min_{(j,k) \in \cS_2} \|f_{j,k}^0\|_2 \ge r_n$ where:

r_n \gtrsim \sqrt{\left(s_1\left(n^{-\frac{2\beta_1}{1 + 2\beta_1}}\log^4{n} + \frac{\log{d}}{n}\right)  + s_2 \left(n^{-\frac{\beta_2}{1 + \beta_2}}\log^4{n} + \frac{\log{d}}{n}\right)\right)} \,.

Note that $r_n$ is precisely the rate of convergence obtained in Corollary <ref>. Assumption <ref> ensures that the signal in the active components is strong enough to be detected with high probability (ensuring the true active set is a subset of the estimated active set). Furthermore, by choosing an optimal penalty level $\lambda_{n, i}$ we control the number of false positives.
Finally, in step 4 of Algorithm <ref>, we optimized squared error loss on the selected subset to produce the final estimate. The following theorem establishes rate of convergence of the estimator $\hat f_{\text{final}}$ obtained via Algorithm <ref>:

Assume that the restricted strong convexity assumption (Assumption <ref>) holds with high probability on $\{X_1, \dots, X_n\}$.
Then, under Assumption <ref>, the estimator $\hat f^{\text{final}}$ satisfies, upto log-factors,

\|\hat f^{\text{final}} - f_0\|_2^2 = O_p\left(s_1\lambda_{n, 1}^2 +  s_2\lambda_{n, 2}^2\right)\,.

where $\lambda_{n, 1}$ and $\lambda_{n, 2}$ are same as in Theorem <ref>.

The proof of the Theorem will be deferred to Appendix <ref>.

Under Assumption <ref> and $\|f^0_j\|_\infty \le b$ for some constant $b$ for all $1 \le j \le p$, with high probability, we have:

\frac{1}{8}\|\hat f^{\text{init}} - f^0\|_n^2 + \frac{\lambda_n}{2} \sum_{j \in \cS^c} \|\hat f^{\text{init}}_j \|_{n} \le C s\lambda_n^2 \,,

where the optimal choice of $\lambda_n \sim n^{-\frac{\beta}{2\beta + 1}} + \sqrt{\frac{\log{p}}{n}}$  .

The detailed proof of the theorem can be found in our main draft. However, here we point out two key ideas which are crucial for the proof of our next theorem:

§.§ Minimax lower bound

In this section, we establish the minimax lower bound on the rate of convergence of the high dimensional two-way interaction model. As mentioned previously, our proofs can be extended to a general $k$-way interaction model. Our main theorem is as follows:

Suppose that we have $n$ observations from model (<ref>). Then under Assumptions <ref>, <ref> we have:

\inf_{\hat f} \sup_{\substack{f \in \mathcal{F}_{\text{sp}} \\ X \sim P_X}} \|\hat f - f\|_2^2 \gtrsim \left(s_1\left(n^{-\frac{2\beta_1}{2\beta_1 + 1}} \vee \frac{\log{\left(d/s_1\right)}}{n} \right) + s_2\left(n^{-\frac{2\beta_2}{2\beta_2 + 2}} \vee \frac{\log{\left(d^2/s_2 \right)}}{n}  \right) \right)

It is immediate from Theorem <ref> and Theorem <ref> that our estimator of the mean function based on neural networks is indeed minimax optimal up to log factors.  Although the proof techniques of Theorem <ref> and Theorem <ref> are similar, there are some important changes in the construction for the alternatives since  we need to enforce sparsity here. We would also like to point out that our proof is different from that of [Raskutti et al., 2012] since we do not assume that the components of the mean function belong to RKHS. Here we sketch the main idea of the proof briefly: each univariate component function is assumed to be $\beta_1$-smooth and each bivariate component function is assumed to be $\beta_2$-smooth. We first construct a collection of alternatives for these components and then we take a sparse combination of them. Finally, constructing alternatives along with a proper choice of the relevant hyper-parameters yield the lower bound. The proof is deferred to the Appendix.

§ CONCLUSION

Deep neural networks have achieved tremendous success in nonparametric function estimation. Yet due to the intrinsic difficulty of “curse-of-dimensionality”, they can only handle nonparametric functions of low-dimension, without excessive restrictions on the function classes. This calls for low-dimensional nonparametric interaction models. At the same time, modern big data applications often involve nonparametric regression with a large number of predictors.  Yet, most statistical theories on neural networks focus only on finite-dimensional regression.  These give rise to the imminent need for the study of nonparametric interaction models in diverging dimensions and understanding the impact of dimensionality in such structured nonparametric models.

This paper contributes critically to understanding the performance of neural networks in low-order interaction models in diverging dimensions.  An important conclusion of our study is that estimated components should have low biases in order to avoid unnecessary error accumulations and this is achieved by our newly introduced debiasing techniques.
For slowly diverging dimensional problems, no additional regularization is needed.  With proper debiasing,direct least-squares estimation on the structured neural networks is shown to achieve a rate of convergence that matches with a newly established minimax low bound, namely, it is optimal.  In a high-dimensional setting, sparsity assumption on the interaction terms is necessary.  We appeal to the penalized least-squares estimation and screening techniques (for random design) and show that the resulting procedure is minimax optimal by establishing a matching lower bound.  Our results provide a comprehensive view of the performance of neural networks for the structured nonparametric model in diverging dimensions, which are critical to modern data science and necessary for neural networks to succeed.

This paper is supported by ONR  N00014-22-1-2340 and the NSF grants DMS-2052926, DMS-2053832, DMS-2210833.

[Antoniadis and Fan, 2001]
Antoniadis, A. and Fan, J. (2001).
Regularization of wavelet approximations.
Journal of the American Statistical Association 96

[Barron, 1993]
Barron, A. R. (1993).
Universal approximation bounds for superpositions of a sigmoidal
IEEE Transactions on Information theory 39

[Bartlett et al., 2019]
Bartlett, P. L., Harvey, N., Liaw, C. and
Mehrabian, A. (2019).
Nearly-tight vc-dimension and pseudodimension bounds for piecewise
linear neural networks.
The Journal of Machine Learning Research 20

[Bauer and Kohler, 2019]
Bauer, B. and Kohler, M. (2019).
On deep learning as a remedy for the curse of dimensionality in
nonparametric regression.
The Annals of Statistics 47 2261–2285.

[Belloni and Chernozhukov, 2013]
Belloni, A. and Chernozhukov, V. (2013).
Least squares after model selection in high-dimensional sparse
Bernoulli 19 521–547.

[Bickel et al., 2009]
Bickel, P. J., Ritov, Y. and Tsybakov, A. B. (2009).
Simultaneous analysis of lasso and dantzig selector .

[Bickel et al., 2010]
Bickel, P. J., Ritov, Y., Tsybakov, A. B.
et al. (2010).
Hierarchical selection of variables in sparse high-dimensional
IMS Collections 6 28.

[Bien et al., 2013]
Bien, J., Taylor, J. and Tibshirani, R. (2013).
A lasso for hierarchical interactions.
The Annals of statistics 41 1111.

[Breiman, 2001]
Breiman, L. (2001).
Random forests.
Machine learning 45 5–32.

[Buja et al., 1989]
Buja, A., Hastie, T. and Tibshirani, R. (1989).
Linear smoothers and additive models.
The Annals of Statistics  453–510.

[Candès and Tao, 2007]
Candès, E. and Tao, T. (2007).
The dantzig selector: Statistical estimation when p is much larger
than n.
The Annals of Statistics 35 2313–2351.

[Chernozhukov et al., 2014]
Chernozhukov, V., Chetverikov, D. and Kato, K.
Gaussian approximation of suprema of empirical processes.
The Annals of Statistics 42 1564–1597.

[Fan and Gijbels, 1996]
Fan, J. and Gijbels, I. (1996).
Local polynomial modelling and its applications.
Chapman & Hall.

[Fan and Gu, 2022]
Fan, J. and Gu, Y. (2022).
Factor augmented sparse throughput deep relu neural networks for high
dimensional regression.
arXiv preprint arXiv:2210.02002 .

[Fan et al., 2022]
Fan, J., Gu, Y. and Zhou, W.-X. (2022).
How do noise tails impact on deep relu networks?
arXiv preprint arXiv:2203.10418 .

[Fan et al., 1998]
Fan, J., Härdle, W. and Mammen, E. (1998).
Direct estimation of low-dimensional components in additive models.
The Annals of Statistics 26 943–971.

[Fan et al., 2020]
Fan, J., Li, R., Zhang, C.-H. and Zou, H.
Statistical foundations of data science.
Chapman and Hall/CRC.

[Fan and Lv, 2008]
Fan, J. and Lv, J. (2008).
Sure independence screening for ultrahigh dimensional feature space.
Journal of the Royal Statistical Society: Series B
(Statistical Methodology) 70 849–911.

[Fan and Song, 2010]
Fan, J. and Song, R. (2010).
Sure independence screening in generalized linear models with
The Annals of Statistics 38 3567–3604.

[Freund et al., 1996]
Freund, Y., Schapire, R. E. et al. (1996).
Experiments with a new boosting algorithm.
In icml, vol. 96. Citeseer.

[Friedman and Stuetzle, 1981]
Friedman, J. H. and Stuetzle, W. (1981).
Projection pursuit regression.
Journal of the American statistical Association 76

[Goodfellow et al., 2016]
Goodfellow, I., Bengio, Y. and Courville, A. (2016).
Deep learning.
MIT press.

[Han and Wellner, 2019]
Han, Q. and Wellner, J. A. (2019).
Convergence rates of least squares regression estimators with
heavy-tailed errors.
The Annals of Statistics 47 2286–2319.

[Hao and Zhang, 2014]
Hao, N. and Zhang, H. H. (2014).
Interaction screening for ultrahigh-dimensional data.
Journal of the American Statistical Association 109

[Hardle et al., 1993]
Hardle, W., Hall, P. and Ichimura, H. (1993).
Optimal smoothing in single-index models.
The Annals of Statistics 21 157–178.

[Hornik et al., 1989]
Hornik, K., Stinchcombe, M. and White, H. (1989).
Multilayer feedforward networks are universal approximators.
Neural networks 2 359–366.

[Horowitz and Mammen, 2007]
Horowitz, J. L. and Mammen, E. (2007).
Rate-optimal estimation for a general class of nonparametric
regression models with unknown link functions.
The Annals of Statistics 35 2589–2619.

[Kohler and Krzyżak, 2005]
Kohler, M. and Krzyżak, A. (2005).
Adaptive regression estimation with multilayer feedforward neural
Nonparametric Statistics 17 891–913.

[Kohler and Krzyżak, 2016]
Kohler, M. and Krzyżak, A. (2016).
Nonparametric regression based on hierarchical interaction models.
IEEE Transactions on Information Theory 63

[Kohler and Langer, 2021]
Kohler, M. and Langer, S. (2021).
On the rate of convergence of fully connected deep neural network
regression estimates.
The Annals of Statistics 49 2231–2249.

[Koltchinskii and Yuan, 2010]
Koltchinskii, V. and Yuan, M. (2010).
Sparsity in multiple kernel learning.
The Annals of Statistics 38 3660–3695.

[Kosorok, 2008]
Kosorok, M. R. (2008).
Introduction to empirical processes and semiparametric

[Kpotufe and Dasgupta, 2012]
Kpotufe, S. and Dasgupta, S. (2012).
A tree-based regressor that adapts to intrinsic dimension.
Journal of Computer and System Sciences 78

[Krizhevsky et al., 2009]
Krizhevsky, A., Hinton, G. et al. (2009).
Learning multiple layers of features from tiny images .

[LeCun, 1998]
LeCun, Y. (1998).
The mnist database of handwritten digits.
http://yann. lecun. com/exdb/mnist/ .

[Lin and Zhang, 2006]
Lin, Y. and Zhang, H. H. (2006).
Component selection and smoothing in multivariate nonparametric
The Annals of Statistics 34 2272–2297.

[Lu et al., 2021]
Lu, J., Shen, Z., Yang, H. and Zhang, S.
Deep network approximation for smooth functions.
SIAM Journal on Mathematical Analysis 53

[Mhaskar, 1996]
Mhaskar, H. N. (1996).
Neural networks for optimal approximation of smooth and analytic
Neural computation 8 164–177.

[Mukherjee et al., 2021]
Mukherjee, D., Banerjee, M. and Ritov, Y. (2021).
Optimal linear discriminators for the discrete choice model in
growing dimensions.
The Annals of Statistics 49 3324–3357.

[Novembre et al., 2008]
Novembre, J., Johnson, T., Bryc, K.,
Kutalik, Z., Boyko, A. R., Auton, A.,
Indap, A., King, K. S., Bergmann, S.,
Nelson, M. R. et al. (2008).
Genes mirror geography within europe.
Nature 456 98–101.

[Raskutti et al., 2012]
Raskutti, G., J Wainwright, M. and Yu, B. (2012).
Minimax-optimal rates for sparse additive models over kernel classes
via convex programming.
Journal of machine learning research 13.

[Ravikumar et al., 2009]
Ravikumar, P., Lafferty, J., Liu, H. and
Wasserman, L. (2009).
Sparse additive models.
Journal of the Royal Statistical Society: Series B
(Statistical Methodology) 71 1009–1030.

[Schmidt-Hieber, 2020]
Schmidt-Hieber, J. (2020).
Nonparametric regression using deep neural networks with relu
activation function.
The Annals of Statistics 48 1875–1897.

[Stone, 1982]
Stone, C. J. (1982).
Optimal global rates of convergence for nonparametric regression.
The Annals of statistics  1040–1053.

[Stone, 1985]
Stone, C. J. (1985).
Additive regression and other nonparametric models.
The Annals of Statistics 13 689–705.

[Stone, 1994]
Stone, C. J. (1994).
The use of polynomial splines and their tensor products in
multivariate function estimation.
The Annals of Statistics 22 118–171.

[Tan and Zhang, 2019]
Tan, Z. and Zhang, C.-H. (2019).
Doubly penalized estimation in additive regression with
high-dimensional data.
The Annals of Statistics 47 2567–2600.

[Tsybakov, 2004]
Tsybakov, A. B. (2004).
Introduction to nonparametric estimation, 2009.
URL https://doi. org/10.1007/b13794. Revised and extended
from the 9.

[Vaart and Wellner, 1997]
Vaart, A. v. d. and Wellner, J. A. (1997).
Weak convergence and empirical processes with applications to
Journal of the Royal Statistical Society-Series A Statistics
in Society 160 596–608.

[Van Der Vaart and Wellner, 1996]
Van Der Vaart, A. W. and Wellner, J. A. (1996).
Weak convergence.
In Weak convergence and empirical processes. Springer,

[Yang and Dunson, 2016]
Yang, Y. and Dunson, D. B. (2016).
Bayesian manifold regression.
The Annals of Statistics 44 876–905.

[Yang and Tokdar, 2015]
Yang, Y. and Tokdar, S. T. (2015).
Minimax-optimal nonparametric regression in high dimensions.
The Annals of Statistics 43 652–674.

[Yarotsky, 2017]
Yarotsky, D. (2017).
Error bounds for approximations with deep relu networks.
Neural Networks 94 103–114.

[Yuan and Lin, 2006]
Yuan, M. and Lin, Y. (2006).
Model selection and estimation in regression with grouped variables.
Journal of the Royal Statistical Society: Series B
(Statistical Methodology) 68 49–67.

[Yuan and Zhou, 2016]
Yuan, M. and Zhou, D.-X. (2016).
Minimax optimal rates of estimation in high dimensional additive
The Annals of Statistics 44 2564–2593.

[Zhang et al., 2014]
Zhang, Y., Wainwright, M. J. and Jordan, M. I.
Lower bounds on the performance of polynomial-time algorithms for
sparse linear regression.
In Conference on Learning Theory. PMLR.

[Zhao et al., 2009]
Zhao, P., Rocha, G. and Yu, B. (2009).
The composite absolute penalties family for grouped and hierarchical
variable selection .

§ PROOFS

§.§ Proof of Theorem <ref>

The proof requires a bias correction step which we describe here in detail. The key idea is to approximate each of the component functions $\{f_j\}, \{f_{jk}\}$ separately by a small neural network and then add them up to approximate $f_0$. Addition of the outputs of two neural network can be performed to another (enlarged) neural network, as depicted in Figure <ref>. The key idea is: if $g_1$ and $g_2$ are two NN, then consider another NN $g$, which operates as follows:

* It takes the input $x$, passes it through, $g_1, g_2, -g_1, -g_2$,

* It then gives output $g(x) = \sigma(g_1(x)) - \sigma(-g_1(x)) + \sigma(g_2(x)) - \sigma(-g_2(x)) = g_1(x) + g_2(x)$.

Therefore if $g_1 \in \cF_{NN}(N_1, L_1), g_2 \in \cF_{NN}(N_2, L_2)$, then $g \in \cF_{NN}(2(N_1 + N_2), (L_1 \vee L_2) + 1)$. Similar analysis also establishes that one can perform sum of $I$ may DNNs via another enlarged DNN; in particular:

\{g_i \in \cF_{NN}(N_i, L_i)\}_{1 \le i \le I}
\implies \sum_{i = 1}^I g_i \in \cF_{N, L} \ \ \ N = 2\sum_i N_i, \ \ L = \vee_i L_i + 1 \,.

Now, we cite a Lemma which will be crucial in proving the approximation error.

Given a smooth function $f \in C^s([0, 1]^d)$ with $s \in \bbN_+$, for any $N, L \in \bbN_+$, there exists a function $\phi$ implemented by a ReLU FNN with width $C_1(N + 2)\log_2(8N)$ and depth $C_2(L + 2) \log_2(4L) + 2d$ such that:

\|\phi - f\|_{L_\infty([0,1]^d)} \le C_3\|f\|_{C^s([0,1]^d)}(NL)^{-2s/d}
where $C_1 = 17s^{d+1}3^dd$, $C_2 = 18s^2$, and $C_3 = 85(s + 1)^d8^s$.   Moreover,

\begin{equation}
\|\phi-f\|_{\mathcal{W}^{1,p}([0,1]^d)} \lesssim (s+1)^d8^s (NL)^{-2(s-1)/d}
\end{equation}

By Lemma <ref>, for all $1 \le j, k \le d$ there exits $g_j, g_{jk} \in \cF_{NN}(C_1N\log{N}, C_2L\log{L})$ such that:

* $\|g_j - f_j\|_\infty \le C_3 \|f_j\|_{C^\beta([0, 1])} (NL)^{-2\beta}$

* $\|g_{jk} - f_{jk}\|_\infty \le C_3 \|f_{jk}\|_{C^\beta([0, 1]^2)} (NL)^{-\beta}$  .

Define the Holder norm earlier. So far we haven't impose the constraint that the integral of $f_j$ is $0$ and the marginal of $f_{jk}$ is $0$. Next we show that, there exists neural network $\tilde g_j$ and $\tilde g_{jk}$ such that they have same approximation error as of $g_j$ and $g_{jk}$ (upto a constant multiple) and they satisfy the constraint that $\tilde g_j$ integrates to $0$ and $\tilde g_{jk}$ has $0$ marginals.

Imposing the first constraint is easy: set $\tilde g_j \leftarrow (g_j - \int g_j)$. Note that $\tilde g_j$ is also a neural network with the same architecture as $g_j$ as subtracting a constant amounts to changing the bias of the last layer. We can bound the approximation error of $f_j$ by $\tilde g_j$ as follows:

\|f_j - \tilde g_j\|_\infty = \left\|f_j - \int f_j - g_j + \int g_j\right\|_\infty \le 2\|f_j - g_j\|_\infty \le 2C_3 \|f_j\|_{C^\beta([0, 1])} (NL)^{-2\beta} \,.

Imposing the second constraint, i.e. marginal of $f_{jk}$ is $0$ is more tricky as we need the marginal function of the neural network approximator to be $0$.

By Lemma <ref>, for any $i,j$ with $i \neq j, i,j \in [d]$, we can obtain RELU FNNs with with $O(N\log N)$ and depth $O(L \log L)$ such that

* $\|f_{ij}-\widetilde{g}_{ij}\|_{\infty} \lesssim (NL)^{-\alpha}$,

* $\|\int \widetilde{g}_{ij}(.,y)dy- \int f(.,y) dy\|_{\infty} \lesssim (NL)^{-2\alpha}$,

* $\|\int \widetilde{g}_{ij}(x,.)dx- \int f(x,.) dx\|_{\infty} \lesssim (NL)^{-2\alpha}$, and

* $\iint f_{ij} = \iint g_{ij}$.

Finally, define $\widetilde{f}=\sum_{i\neq j} \widetilde{g}_{ij}$. Then,

\begin{align*}
&\bbE\left(\sum_{i \neq j}f_{ij}-\widetilde{f}\right)^2 = \bbE \left(\sum_{i \neq j}f_{ij} -\sum_{i \neq j} \widetilde{g}_{ij}\right)^2 \leq p_{\max} \iint \left(\sum_{i \neq j}f_{ij} -\sum_{i \neq j} \widetilde{g}_{ij}\right)^2 d\mathbf{x}\\
&=p_{\max}\left(\sum_{i\neq j} \iint (f_{ij}-\widetilde{g}_{ij})^2 dx_idx_j+ \sum_{i \neq j\neq k}\int \left(\int f_{ij}-\widetilde{g}_{ij}dx_j\right)\left(\int f_{ik}-\widetilde{g}_{ik}dx_k\right)dx_i\right) \\
&\le p_{max} \Big(d^2 (NL)^{-2\alpha} +d^3 (NL)^{-4\alpha}\Big),
\end{align*}

by our construction of $\widetilde{g}_{ij}$. Need to assume cov term smaller

§.§ Proof of Theorem <ref>

We want to invoke <cit.>. For any $\phi = \phi_1 + \phi_2$ (where $\phi_1 \in \cF_{NN}^1, \phi_2 \in \cF^2_{NN}$),

define the population risk function $R(\phi)$ (resp. empirical risk $\hat R_n(\phi)$) as

R(\phi) = \bbE[(Y - \phi(X))^2] \ \ \ \left(\text{resp. } \hat R(\phi) = \frac1n \sum_i (Y_i - \phi(X_i))^2\right)

Since that $\eps$ is independent of $X$ with mean $0$, it is immediate that:

R(\phi) - R(f_0) = \bbE[(f(X) - f_0(X))^2] \triangleq d^2(\phi, f_0) \,.

Define $\cF_n := \{\phi = \phi_1 + \phi_2: \phi_1 \in \cF^1_{NN}, \phi_2 \in \cF^2_{NN}\}$ (i.e. $\cF_n = \cF_{NN}^1 + \cF_{NN}^2$) and $\omega_n = d(\phi^\star, f_0)$, i.e. the approximation error of the class of neural network. It is immediate that:

\textstyle
\cF_n \subseteq \cF_{NN}\left(d, c_1 \left(dN_1\log{N_1} + \binom{d}{2}N_2\log{N_2}\right), c_2, W_n, 1\right)

W_n = c_3\left(d(N_1\log{N_1})^2 + \binom{d}{2}\left(N_2\log{N_2}\right)^2\right) \,,

for some constant $c_1, c_2, c_3$ independent of $(n, d)$. Therefore, VC dimension of $\cF_n$ is less than or equal to $V_n \triangleq c_4 W_n \log{W_n}$ for some $c_4>0$. We use this in the subsequent analysis. Let $\zeta_n(\delta)$ be such that:

\begin{equation}
\label{eq:max_ineq_1}
\bbE\left[\sup_{\phi \in \cF_n: d(\phi, f_0) \le \delta} \left|\left(R_n(\phi) - R(\phi)\right) - \left(R_n(f_0) - R(f_0)\right)\right|\right] \lesssim \frac{\zeta_n(\delta)}{\sqrt{n}} \,.
\end{equation}

This function $\zeta_n$ is called the modulus of continuity, which is used to bound the fluctuation of the empirical process around the population process in a neighborhood of the true function. $\zeta_n(\delta)$ intrinsically depends on the complexity of the underlying function class. We quantify this dependency through the following maximal inequality (e.g. Theorem 5.2 of [Chernozhukov et al., 2014]), which we state here for the convenience of the readers:

Define the operator $\bbG_n = \sqrt{n}(\bbP_n - P)$. Consider a collection of functions $\cF$ with envelope function $F$. Assume that $F \in L_2(P)$. Then:

\bbE\left[\sup_{f \in \cF} \left| \bbG_n(f)\right|\right] \lesssim J(r, \cF, F)\|F\|_2 + \frac{\|M\|_2 J^2(r, \cF, F)}{r^2 \sqrt{n}}

for $r = \sigma/\|F\|_2$,  where

\begin{align*}
J(\tau, F, \cF) & = \int_0^\tau \sup_Q \sqrt{1 + \log{N\left(\eps \|F\|_2, \cF, L_2(Q)\right)}} d \eps \\
M & = \max_{1 \le i \le n} F(X_i), \qquad  \sigma^2 = \sup_{f \in \cF} Pf^2   \,.
\end{align*}

We use Lemma <ref> to find the function $\zeta_n(\delta)$.

For any $\phi$, define a function $g \equiv g(\phi)$ on the space of $(X, \eps)$ as:

\begin{align}
g(X, \eps) & = (f_0(X) - \phi(X) + \eps)^2 - \eps^2 \notag \\
\label{eq:g_val} & = (\phi(X) - f_0(X))^2 + 2(\phi(X) - f(X))\eps.
\end{align}

Define the collection of $\cG$ as $\cG = \left\{g(\phi) : \phi \in \cF_{n} \right\}$. This definition implies:

\bbE\left[\sup_{\phi: d(\phi, f_0) \le \delta} \left|\left(R_n(\phi) - R(\phi)\right) - \left(R_n(f_0) - R(f_0)\right)\right|\right] = \frac{1}{\sqrt{n}}\bbE\left[\sup_{g \in \cG_\delta} \left|\bbG_n(g)\right|\right].

Here we $\cG_\delta$ as the collection of all $g(\phi)$'s, such that $d(\phi, f_0) \le \delta$.
To apply Lemma <ref>, we need to i) quantify $\sigma^2$, ii) find an envelope of $\cG_\delta$ and iii) bound the logarithm of covering number. For i) observe that:

\begin{align*}
\sup_{g \in \cG: d(f, f_0) \le \delta} \bbE[g(X)^2] \le \delta^2(2 + 8\sigma_\eps^2)  \triangleq C_{\sigma_\eps}\delta^2 \,.
\end{align*}

Next, to construct an envelope, note that (<ref>) implies:

|g(X, \eps)| \le 9B^2 + 6B\eps \,.

for all $(X, \eps)$. This follows from the fact that $\|\phi\|_\infty \le 2B$ and $\|f_0\|_\infty \le B$. Therefore, we can take the envelope $G$ as $G(X, \eps) = 9B^2 + 6B\eps$.

To bound the logarithm of the covering number of $\cG_\delta$,

we use Lemma 9.9 of [Kosorok, 2008]. Note that, each $g=g_1 + g_2 + h g_3$ where $g_1 = (\phi - f_0)^2_+$, $g_2 = (\phi - f_0)^2_-$, $g_3 = 2(\phi - f_0)$ and $h(x, \eps) = \eps$. We calculate the VC dimensions of each class of functions:

* The VC dimension of $\cG_1 = \left\{g_1(\phi): \phi \in \cF_{NN}\right\}$ is $\le V_n$. To see this, write $g_1(x, \eps) = m_1 \circ (\phi - f_0) \circ m_2(x, \eps)$ where $m_2(x, \eps) = x$ and $m_1(x) = x_+^2$. As precomposing by a fixed function $\phi$, subtracting a fixed function $f_0$ and post-composing by a monotone function does not change the VC dimension, the claim follows. Similar argument establishes that VC dimension of $\cG_2 = \left\{g_2(f): f \in \cF\right\}$ is $\le V_n$.

* VC dimension of $\cG_3$, the collection of $g_3h$ is $\lesssim V_n$. To see this, note that VC dimension of $g_3$ is $\le V_n$. As before, $g_3 = 2 (\phi - f_0) \circ \psi$ as we have already mentioned that pre-composing by $\psi$, subtracting $f_0$ and multiplying by by a fixed function (here $2h$) does not change the VC dimension.

Therefore, we have established that VC-dim of $\cG_i$ is $\lesssim V_n$. This implies, via Hausslers's bound <cit.>):

\sup_Q \log{N\left(\eps \|G_i\|_2, \cG_i, L_2(Q)\right)} \lesssim
\log{KV_n} + 2V_n\log{\left(\frac{2\sqrt{2e}}{\eps}\right)} \,.

As $\cG_\delta \subseteq \cG_1 + \cG_2 + \cG_3$, we have:

\begin{align*}
\log{N\left(\eps \|G\|_2, \cG_\delta, L_2(Q)\right)} & \lesssim \sum_{j = 1}^3 \log{N\left(\eps \|G_i\|_2, \cG_i, L_2(Q)\right)} \\
& \le 3\log{KV_n} + 6V_n\log{\left(\frac{2\sqrt{2e}}{\eps}\right)} \\
& \le 9V_n\log{\left(\frac{2\sqrt{2e}}{\eps}\right)} \,,
\end{align*}

where the last inequality holds as soon as $\log{KV_n} \le V_n$ and $\eps \le 2\sqrt{2/e}$. Therefore, we can use the last inequality for all small $\eps$ and for all large enough $V_n$. Taking $r = \sqrt{C_{\sigma_\eps}}\delta/\|G\|_2 $ we have:

\begin{align*}
J(r, \cG_\delta, G) & \le \int_0^r \sqrt{1 + 9V_n\log{\left(\frac{2\sqrt{2e}}{\eps}\right)}} \ d\eps \\
& \le 4\sqrt{V_n}\int_0^r \sqrt{\log{\left(\frac{2\sqrt{2e}}{\eps}\right)}} \ d\eps \\
& = 8\sqrt{2eV_n}\int_0^{r/2\sqrt{2e}} \sqrt{\log{\left(\frac{1}{\eps}\right)}} \ d\eps \\
& \le 4r\sqrt{V_n}\sqrt{\log{\left(\frac{2\sqrt{2e}}{r}\right)}}.
\end{align*}

Here the second inequality holds as long as $V_n\log{(2\sqrt{2}/\eps)} \ge 1/7$, i.e. $\eps < 2\sqrt{2}\exp{\left(-1/(7V_n)\right)}$. This holds for all large $n$ as long as $r < 1$. Hence we conclude via Lemma <ref>:

\begin{align*}
\bbE\left[\sup_{g \in \cG_\delta} \left|\bbG_n(g)\right|\right] & \lesssim  \delta \sqrt{V_n\log{\left(\frac{1}{\delta}\right)}} + \frac{V_n}{\sqrt{n}}\log{\left(\frac{1}{\delta}\right)}\sqrt{\bbE\left[\max_i (1 + 2|\eps_i|)^2\right]}  \\
& \lesssim  \delta \sqrt{V_n\log{\left(\frac{1}{\delta}\right)}} + \frac{V_n \sqrt{\log{n}}}{\sqrt{n}}\log{\left(\frac{1}{\delta}\right)} \triangleq \phi_n(\delta) \,.
\end{align*}

Now we establish the rate of convergence. Set $r_n = \left(\sqrt{\frac{V_n}{n}\log{\frac{n}{V_n} }\sqrt{\log{n}}} + \omega_n\right)^{-1}$, where $\omega_n = d(\phi^\star, f_0)$ is the approximation error defined before. Note that $r_n$ satisfies the equation $\phi_n(r_n^{-1}) \lesssim r_n^{-2}\sqrt{n}$ since

\begin{align*}
r_n^2\phi_n\left(\frac{1}{r_n}\right) & \lesssim r_n\sqrt{V_n \log{r_n}} + r_n^2 \frac{V_n \sqrt{\log{n}}}{\sqrt{n}}\log{r_n} \\
& \le \sqrt{n}\left[\left(\frac12\frac{\log{\frac{n}{V_n} - \log{\log{\frac{n}{V_n}}} - \frac12 \log{\log{n}}}}{\log{\frac{n}{V_n} }\sqrt{\log{n}}}\right)^{1/2} + \left(\frac12\frac{\log{\frac{n}{V_n} - \log{\log{\frac{n}{V_n}}} - \frac12 \log{\log{n}}}}{\log{\frac{n}{V_n} }}\right)\right] \\
& \lesssim \sqrt{n}
\end{align*}

We now use shelling argument to establish the rate of convergence. Fix $t > 1$, and define $A_j := \{\cF_{NN} : 2^{j-1}t \le r_n d(f_0, \phi) \le 2^jt\}$.  It follows from the definition of ERM, we have

\begin{align}
& \bbP\left(r_n d(\hat \phi, f_0) \ge t\right) \notag = \bbP\left(\sup_{\substack{\phi \in \cF_{NN} \\ r_n d(f_0, \phi) \ge t}} \hat R_n(\phi^\star) - \hat R_n(\phi) \ge 0 \right) \notag \\
& \le \sum_{j=1}^\infty \bbP\left(\sup_{\phi \in A_j} \hat R_n(\phi^\star) - \hat R_n(\phi) \ge 0 \right) \notag \\
%& \triangleq \sum_{j=1}^\infty \bbP\left(\sup_{\phi \in A_j} \hat R_n(\phi^\star) - \hat R_n(\phi) \ge 0 \right)  \notag \\
& \le \sum_{j=1}^\infty \bbP\left(\sup_{\phi \in A_j} \left(R_n(\phi^\star) - R(\phi^\star)\right) - \left(R_n(\phi) - R(\phi)\right) \ge \inf_{\phi \in A_j} R(\phi) - R(\phi^\star) \right) \notag \\
\label{eq:rate_1} & \le \sum_{j=1}^\infty \frac{\bbE\left[\sup_{\phi \in A_j} \left|\left(R_n(\phi^\star) - R(\phi^\star)\right) - \left(R_n(\phi) - R(\phi)\right)\right|\right]}{\inf_{\phi \in A_j} R(\phi) - R(\phi^\star) }
\end{align}

We bound the numerator and denominator separately. For the denominator note that:

\begin{align*}
\inf_{f \in A_j} R(\phi) - R(\phi^\star)  & \ge \inf_{f: d_{f, f_0} \ge 2^{j-1}tr_n^{-1}} R(\phi) - R(f_0) + R(f_0) - R(\phi^\star) \\
& \ge  2^{2j-2}t^2r_n^{-2} - \omega^2_n
\end{align*}

For the numerator:

\begin{align*}
& \bbE\left[\sup_{\phi \in A_j} \left|\left(R_n(\phi^\star) - R(\phi^\star)\right) - \left(R_n(\phi) - R(\phi)\right)\right|\right] \\
& \le \bbE\left[\sup_{f: d(f, f_0) \le 2^j tr_n^{-1}} \left|\left(R_n(\phi^\star) - R(\phi^\star)\right) - \left(R_n(\phi) - R(\phi)\right)\right|\right] \lesssim \frac{\zeta_n\left(2^j tr_n^{-1}\right)}{\sqrt{n}} \,.
\end{align*}

Here, the second last equation holds because $d(\phi^\star, f_0) = \omega_n \le 2^j tr_n^{-1}$ as $r_n^{-1} \ge \omega_n$. Putting this bound in equation (<ref>) we have:

\begin{align*}
\bbP\left(r_n d(\hat f, f_0) \ge t\right) & \le \sum_{j=1}^\infty \frac{\phi_n\left(2^j tr_n^{-1}\right)}{\sqrt{n}\left(2^{2j-2}t^2r_n^{-2} - \omega_n^2\right) } \\
& \le \sum_{j=1}^\infty \frac{2^j t\phi_n\left(r_n^{-1}\right)}{\sqrt{n}\left(2^{2j-2}t^2r_n^{-2} - \omega^2_n \right) } =  \sum_{j=1}^\infty \frac{2^j t r_n^2 \phi_n(r_n^{-1})n^{-1/2}}{2^{2j-2}t^2 - \omega^2_nr_n^2 } \\
& \lesssim \sum_{j=1}^\infty \frac{2^j t}{2^{2j-2}t^2 - 1 } \hspace{0.2in} [\text{since } r_n^2 \phi_n(r_n^{-1}) \lesssim \sqrt{n}, \ \ \omega_n^2r_n^2 \le 1]\\
& = \frac{1}{t} \sum_{j=1}^\infty \frac{2^j}{2^{2j-2} - \frac{1}{t^2} } \le \frac{c}{t}
\end{align*}

This proves that:

r_n \ d\left(\hat \phi, f_0\right) = O_p(1) \,.

\|\hat f - f_0\|_{L_2(P_X)}^2 = O_p\left(\omega_n^2 + \frac{V_n}{n}\log^{3/2}{n}\right) \,.

Now to balance $\omega_n$ and $V_n$, we choose $N_1 = \lfloor n^{1/2(2\beta_1 + 1)}\rfloor$ and $N_2 = \lfloor n^{1/2(\beta_2 + 1)}\rfloor$. This implies:

\begin{align*}
\omega_n^2 & \lesssim dn^{-\frac{2\beta_1}{2\beta_1 + 1}} + \binom{d}{2}n^{-\frac{\beta_2}{\beta_2 + 1}} \\
\frac{V_n}{n} & \lesssim \frac{W}{n}\log{W} \lesssim \left(dn^{-\frac{2\beta_1}{2\beta_1 + 1}} + \binom{d}{2}n^{-\frac{\beta_2}{\beta_2 + 1}}\right)\log^3{n},
\end{align*}

which yields

\|\hat f - f_0\|_{L_2(P_X)}^2 = O_p\left(\left(dn^{-\frac{2\beta_1}{2\beta_1 + 1}} + \binom{d}{2}n^{-\frac{\beta_2}{\beta_2 + 1}}\right)\log^{4.5}{n}\right).

This completes the proof.

§.§ Proof of Theorem <ref>

Here, we establish the lower bound on the following two-way interaction model:

Y_i = \mu + \sum_{j = 1}^d f_j(X_{ij}) + \sum_{k < l} f_{kl}(X_{ik}, X_{il}) + \eps_i \triangleq f_0(X_i) + \eps_i \,.

We show that:

\inf_{\hat f} \sup_{f \in \cF, X \sim P_{X}} \bbE_{P_{X}}[(\hat f(X) - f(X))^2] \ge c\left(dn^{-\frac{2\beta_1}{2\beta_1 + 1}} + \dbinom{d}{2} n^{-\frac{\beta_2}{\beta_2 + 1}}\right) \,.

where $\cF$ is the collection of all additive functions which satisfies: I) The one-dimensional components are $\beta_1$-smooth, $[0, 1]$ and integrate to $0$ and II) the two-dimensional components are $\beta_2$-smooth, supported on $[0, 1]^2$ having marginals $0$. Our proof is based on the techniques introduced in <cit.>. Throughout the proof, assume $X_{ij}\overset{\mathrm{iid}}{\sim} Unif(0,1)$, $i\le n$, $j \le d$ and $\varepsilon_i \overset{\mathrm{iid}}{\sim} N(0,1)$, $i
\le n$. The key idea is to use Fano's inequality upon carefully choosing a subset of $\cF$.

Following [Tsybakov, 2004], fix the following notations:

\begin{equation}\label{eq:define_mi}
m_i = \lfloor c_0 n^{\frac{1}{2\beta_i + 1}} \rfloor, \quad h_i = m_i^{-1}, \quad i=1,2.
\end{equation}

We construct the alternatives in two steps, for one-dimensional components and for two-dimensional components.

Alternatives for one-dimensional components:  First divide $[0, 1]$ into $m_1$ grids $[(k-1)/m_1, k/m_1]$, $1 \le k \le m_1$. For each $k$, define $x_k = (k - 0.5)/m_1$. Select a smooth kernel $K \in C^\infty(\reals)$ supported on $(-1/2, 1/2)$, $\int K(u) \ du = 0$ and define $\int K^2(u) \ du \triangleq \|K\|_2^2$. Define the functions $\phi_k(x)$, $k \le m_1$ as:

\begin{equation}\label{eq:define_lower_phi_lin}
\phi_k(x) := Lh_1^{\beta_1}K\left(\frac{x - x_k}{h_1}\right)
\end{equation}

The constant $L$ is chosen so that $\phi_k$ has bounded derivatives and satisfies Assumption <ref>. For each $k$, the function $\phi_k$ is supported on $[(k-1)/m, k/m]$, and the integral of $\phi_k$ is $0$, since

\begin{align*}
\int \phi_k(x) \ dx & = Lh_1^{\beta_1}\int K\left(\frac{x - x_k}{h_1}\right) \ dx \\
&=  Lh_1^{\beta_1}\int_{x_k - \frac{h}{2}}^{x_k + \frac{h}{2}} K\left(\frac{x - x_k}{h_1}\right) \ dx \hspace{0.2in} [\text{Since } K(z) > 0 \text{ when } |z| \ge 1/2]\\
& = Lh_1^{\beta_1 + 1} \int_{-\frac12}^{\frac12} K(z) \ dz = 0 \,.
\end{align*}

For $j \neq k$, $\phi_j$ and $\phi_k$ have disjoint support, yielding the $L_2([0, 1])$ inner product between any $\phi_j$ and $\phi_k$ is $0$.

Set $M_1 := dm_1$ and $\Omega_1 := \{0, 1\}^{M_1}$. For each $\bomega_1 \in \Omega_1$, define its entries by $\{\bomega_{1, jk}\}_{1 \le j \le d, 1 \le k \le m_1}$. Define functions $f^{(1)}_{\omega_1}$ as:

f^{(1)}_{\bomega_1}(X) := \sum_{j=1}^d \sum_{k=1}^m \bomega_{1, jk}\phi_k(X_j) \,,

which will contribute to the one dimensional alternatives.

Alternatives for two-dimensional components:  The construction for two-dimensional alternatives are similar to that for one-dimensional components. As before, divide $[0, 1]^2$ into $m_2^2$ rectangles $[(k-1)/m_2, k/m_2] \times [(l-1)/m_2, l/m_2]$, $1 \le k, l \le m_2$.
Set $x_k := (k - 0.5)/m_2$. Using $K$ same as before, define $\phi_{kl}$, $1 \le k,l \le m_2$ as

\begin{equation}\label{eq:define_lower_phi_quad}
\phi_{kl}(x, y) := Lh_2^{\beta_2}K\left(\frac{x - x_k}{h}\right)K\left(\frac{y - x_l}{h}\right)
\end{equation}

Observe that, each $\phi_{kl}$ is supported on $[(k - 1)/m_2, k/m_2] \times [(l-1)/m_2, l/m_2]$. and consequently $L_2([0,1]^2)$ inner-product of any two $\phi_{k_1, l_1}$ and $\phi_{k_2, l_2}$ is 0 if $(k_1,l_1)\neq (k_2,l_2)$. Also the marginals of $\phi_{kl}$ are $0$ as $\int K((x - x_k)/h) \ dx = 0$.

Set $M_2 := d(d-1)m^2_2/2$ and $\Omega_2 := \{0, 1\}^{M_2}$. For each $\bomega_2 \in \Omega_2$ we write it in a tensor form $\bomega_{2, ijkl}$ where $1 \le i < j \le d$ and $1 \le k \le l \le m_2$ and define $f^{(2)}_{\bomega_2}$ as:

f^{(2)}_{\bomega_2}(X)= \sum_{i < j} \sum_{k, l} \bomega_{2, ijkl}\phi_{kl} (X_i, X_j)\,.

These functions will be the two-dimensional component of our alternatives.

Final step:  We now construct our alternatives by combining the one-dimensional and two-dimensional components as constructed above. We first set the true mean function $f_0 = 0$. Now we choose $S_1 \subset \Omega_1$ and $S_2 \subset \Omega_2$ carefully and then construct our alternatives as:

\cF_S = \left\{f_\bomega = f^{(1)}_{\bomega_1} + f^{(2)}_{\bomega_2}, \ \bomega_1 \in S_1, \bomega_2 \in S_2 \right\}

To choose $S_1$ and $S_2$ we invoke Varshamov-Gilbert theorem, cf. <cit.> which we state here for the convenience of the reader:

For any $M\ge 8$, there exists $S \subseteq \{0, 1\}^M$ with $|S| \ge 2^{M/8}$, such that for any $\bomega \neq \bomega' \in S$, $\rho(\bomega, \bomega') \ge M/8$, where $\rho$ is the Hamming distance.

By Proposition <ref>, we can choose $S_1$ and $S_2$ such that for any $\bomega_1, \bomega_1' \in S_1$ we have $\rho(\bomega_1, \bomega_1') \ge M_1/8$ and for any $\bomega_2, \bomega_2' \in S_2$, $\rho(\bomega_2, \bomega_2') \ge M_2/8$. Observe that for any $\bomega = (\bomega_1, \bomega_2) \neq \bomega' = (\bomega'_1, \bomega_2')$ we have:

\begin{align*}
& d^2\left(f_\bomega, f_{\bomega'}\right) = \int_{[0, 1]^d} (f_\bomega(\bx) - f_{\bomega'}(\bx))^2 \ d\bx  \\
& = \int_{[0, 1]^d} \left(\sum_{j=1}^d \sum_{k=1}^m (\bomega_{1, jk} - \bomega'_{1, jk})\phi_k(X_j) + \sum_{i < j} \sum_{k, l} (\bomega_{2, ijkl} - \bomega'_{2, ijkl})\phi_{kl}(X_i, X_j)\right)^2 \ dx \\
& = \int_{[0, 1]^d} \left(\sum_{j=1}^d \sum_{k=1}^m (\bomega_{1, jk} - \bomega'_{1, jk})\phi_k(X_j) \right)^2 \ d\bx + \left(\sum_{i < j} \sum_{k, l} (\bomega_{2, ijkl} - \bomega'_{2, ijkl})\phi_{kl}(X_i, X_j)\right)^2 \ d\bx \\
& \qquad \qquad + 2\int_{[0, 1]^d} \left(\sum_{j=1}^d \sum_{k=1}^m (\bomega_{1, jk} - \bomega'_{1, jk})\phi_k(X_j) \right)\left(\sum_{i < j} \sum_{k, l} (\bomega_{2, ijkl} - \bomega'_{2, ijkl})\phi_{kl}(X_i, X_j)\right) \ d\bx \\
& \triangleq T_1 + T_2 + 2T_3 \,.
% & = \int_{[0, 1]^d} \left(\sum_{ijkl}\omega^{ij}_{kl}\phi_{kl}(x_i, x_j) - \sum_{ijkl} (\omega')^{ij}_{kl} \phi_{kl}(x_i, x_j)\right)^2 \ dx \\
% & =  \sum_{ijkl}\left(\omega^{ij}_{kl} - (\omega')^{ij}_{kl} \right)^2 \int_{[0, 1]^d} \phi^2_{kl}(x_i, x_j) \ dx  \hspace{0.2in} [\text{Disjoint support of }\phi_{kl}'s]\\
% & = \rho(\omega, \omega') \int_{[0, 1]^2} \phi^2_{kl}(x_i, x_j) \ dx_i dx_j \\
% & = \rho(\omega, \omega')  L^2 h^{4\alpha + 2} \left(\int_{1/2}^{1/2}K^2(x) dx\right)^2 := L^2\|K\|^4 h^{4\beta + 2} \rho(\omega, \omega') \,.
\end{align*}

We now analyze each $T_i$ separately.

\begin{align*}
T_1 &= \int_{[0, 1]^d} \left(\sum_{j=1}^d \sum_{k=1}^m (\bomega_{1, jk} - \bomega'_{1, jk})\phi_k(X_j) \right)^2 \ d\bx \\
& = \sum_{j=1}^d \sum_{k=1}^m (\bomega_{1, jk} - \bomega'_{1, jk})^2\int_0^1 \phi^2_k(X_j) \ dX_j \\
& = L^2 h_1^{2\beta_1}\sum_{j=1}^d \sum_{k=1}^m (\bomega_{1, jk} - \bomega'_{1, jk})^2\int_0^1 K^2\left(\frac{X_j - x_k}{h_1}\right) \ dX_j \\
& = L^2 h_1^{2\beta_1 + 1}\sum_{j=1}^d \sum_{k=1}^m (\bomega_{1, jk} - \bomega'_{1, jk})^2\int_{-\frac12}^{\frac12} K^2\left(z\right) \ dz \\
& = L^2 h_1^{2\beta_1 + 1}\|K\|_2^2 \ \sum_{j=1}^d \sum_{k=1}^m (\bomega_{1, jk} - \bomega'_{1, jk})^2 = L^2 h_1^{2\beta_1 + 1}\|K\|_2^2 \ \rho(\bomega_1, \bomega'_1)
\end{align*}

For $T_2$:

\begin{align*}
T_2 & = \int_{[0, 1]^d}\left(\sum_{i < j} \sum_{k, l} (\bomega_{2, ijkl} - \bomega'_{2, ijkl})\phi_{kl}(X_i, X_j)\right)^2 \ d\bx \\
& = \sum_{i < j} \sum_{k, l} (\bomega_{2, ijkl} - \bomega'_{2, ijkl})^2 \int_{[0, 1]^2} \phi^2_{kl}(X_i, X_j) \ dX_i \ dX_j \\
& = L^2h_2^{2\beta_2}\sum_{i < j} \sum_{k, l} (\bomega_{2, ijkl} - \bomega'_{2, ijkl})^2 \int_{[0, 1]^2} K^2\left(\frac{X_i - x_j}{h_2}\right)K^2\left(\frac{X_j - x_k}{h_2}\right) \ dX_i \ dX_j \\
& =  L^2h_2^{2\beta_2 + 2} \sum_{i < j}\sum_{k, l} (\bomega_{2, ijkl} - \bomega'_{2, ijkl})^2 \left(\int_{-\frac12}^{\frac12} K^2(z) \ dz\right)^2 \\
& = L^2h_2^{2\beta_2 + 2}\|K\|_2^4 \ \sum_{i < j}\sum_{k, l} (\bomega_{2, ijkl} - \bomega'_{2, ijkl})^2 = L^2h_2^{4\beta_2 + 2}\|K\|_2^4 \rho(\bomega_2, \bomega_2') \,.
\end{align*}

Finally we claim that $T_3 = 0$. This follows from the fact that if $i, j, k$ all are distinct, then for any $1 \le v_1 \le m_1$ and $1 \le v_2, v_3 \le m_2$:

\begin{align*}
& \int_{[0, 1]^3}\phi_{v_1}(X_i) \phi_{v_2, v_3}(X_j, X_k) \ dX_i \ dX_j \ dX_k \\
& = \int_{[0, 1]}\phi_{v_1}(X_i) \ dX_i \ \int_{[0, 1]^2}\phi_{v_2, v_3}(X_j, X_k) \ dX_j \ dX_k = 0
\end{align*}

as the integrals of $\phi_{v_1}$ and $\phi_{v_2, v_3}$ are $0$ by construction. When they are not all distinct, suppose $i =j \neq k$. Then:

\begin{align*}
& \int_{[0, 1]^2}\phi_{v_1}(X_i) \phi_{v_2, v_3}(X_i, X_k) \ dX_i \ dX_k \\
& = \int_{[0, 1]}\phi_{v_1}(X_i) \ \left(\int_{[0, 1]}\phi_{v_2, v_3}(X_i, X_k) \ dX_k\right) \ dX_i = 0
\end{align*}

as marginals of $\phi_{v_2, v_3}$ is 0 by our construction. Therefore $T_3 = 0$. Combining the expression of $T_i$'s, we have:

\begin{align}
\label{eq:d_val}
d^2(f_\bomega, f_{\bomega'}) & =  L^2 h_1^{2\beta_1 + 1}\|K\|_2^4 \ \rho(\bomega_1, \bomega_1') +  L^2h_2^{2\beta_2 + 2} \|K\|_2^4\ \rho(\bomega_2, \bomega_2')
\end{align}

where $\rho$ is the Hamming distance.

From Proposition <ref>, we have:

\begin{align}
\label{eq:d_val_lb}
d^2(f_\bomega, f_{\bomega'}) & \ge  L_1^2 h_1^{2\beta_1 + 1}\|K\|_2^2 \ \frac{dm_1}{8} +  L_2^2h_2^{2\beta_2 + 2} \|K\|_2^4\ \frac{d(d-1)m_2^2}{16} \notag \\
& =  \frac{L_1^2 \|K\|_2^2}{8} dm_1^{-2\beta_1} +  \frac{L_2^2 \|K\|_2^4}{8}\dbinom{d}{2}m_2^{-2\beta_2} := 4\delta^2
% & \ge c\left(\frac{L_1^2 \|K\|_2^2}{8} dn^{-\frac{2\beta_1}{2\beta_1 + 1}} +  \frac{L_2^2 \|K\|_2^4}{8}\dbinom{d}{2}n^{-\frac{2\beta_2}{2\beta_2 + 2}}\right)
\,.
\end{align}

Here the second equality follows from the fact that $h_1 = m_1^{-1}$ and $h_2= m_2^{-1}$ and the third inequality follows from the definition of $m_1, m_2$.
Applying Fano's inequality <cit.> on the collection $\cF_S$, we obtain:

\begin{equation}
\label{eq:Fano}
\inf_{\hat f} \sup_{f,P_{X}} \bbE_{P_{X}}[(\hat f(X) - f(X))^2] \ge \delta^2 \left(1 -\frac{\frac{n}{|S|^2} \sum_{\omega \neq \omega' \in S}KL\left(\bbP_\omega | \bbP_{\omega'}\right) + \log{2}}{\log{(|S| - 1)}}\right)
\end{equation}

As the error $\eps_i$'s are normal, we have:

\begin{align*}
KL\left(\bbP_\omega | \bbP_{\omega'}\right) & = \frac{1}{2}\bbE\left[\left(f_\omega(X) - f_{\omega'}(X)\right)^2\right] = \frac12 d^2\left(f_\bomega, f_{\bomega'}\right) \\
& =  L^2 h_1^{2\beta_1 + 1}\|K\|_2^2 \ \rho(\bomega_1, \bomega_1') +  L^2h_2^{2\beta_2 + 2} \|K\|_2^4\ \rho(\bomega_2, \bomega_2') \\
& \le L^2 h_1^{2\beta_1 + 1}\|K\|_2^2 \ M_1+  L^2h_2^{2\beta_2 + 2} \|K\|_2^4\ M_2 \\
& = L^2 h_1^{2\beta_1 + 1}\|K\|_2^2 \ dm_1+  L^2h_2^{2\beta_2 + 2} \|K\|_2^4\ \dbinom{d}{2} m_2^2 \\
& = L^2\|K\|_2^2 dm_1^{-2\beta_1}+  L^2\|K\|_2^4\ \dbinom{d}{2} m_2^{-2\beta_2} \le 32 \delta^2 \,.
\end{align*}

Moreover, $|S| = |S_1 \times S_2| \ge 2^{(M_1 + M_2)/8}$ by Proposition <ref>. Plugging the bounds in (<ref>), we have:

\begin{align}
\label{eq:FANO_2}
\inf_{\hat f} \sup_{f,P_{X}} \bbE_{P_{X}}[(\hat f(X) - f(X))^2] \ge \delta^2 \left(1 -9\frac{32n\delta^2 + \log{2}}{M_1 + M_2}\right)
\end{align}

Now from the definition of $m_1, m_2$ we have:

\begin{align*}
n\delta^2 & = n\left(\frac{L_1^2 \|K\|_2^2}{32} dm_1^{-2\beta_1} +  \frac{L_2^2 \|K\|_2^4}{32}\dbinom{d}{2}m_2^{-2\beta_2}\right) \\
& \le Cn \left(\frac{L_1^2 \|K\|_2^2}{32} dn^{-\frac{2\beta_1}{2\beta_1 + 1}} +  \frac{L_2^2 \|K\|_2^4}{32}\dbinom{d}{2}n^{-\frac{2\beta_2}{2\beta_2 + 2}}\right)  \\
& = C\left(\frac{L_1^2 \|K\|_2^2}{32} \vee \frac{L_2^2 \|K\|_2^4}{32}\right)\left( dn^{\frac{1}{2\beta_1 + 1}} +  \dbinom{d}{2}n^{\frac{2}{2\beta_2 + 2}}\right) \\
& \triangleq C\left(\frac{L_1^2 \|K\|_2^2}{32} \vee \frac{L_2^2 \|K\|_2^4}{32}\right)\psi_n \,.
\end{align*}

On the other hand:

\begin{align*}
M_1 + M_2 & = dm_1 + \dbinom{d}{2}m_2^2 \ge c\left(dn^{\frac{1}{2\beta_1 + 1}} +  \dbinom{d}{2}n^{\frac{2}{2\beta_2 + 2}}\right) = c\psi_n \,.
\end{align*}

Using these bounds in (<ref>) yields:

\begin{align*}
\inf_{\hat f} \sup_{f,P_{X}} \bbE_{P_{X}}[(\hat f(X) - f(X))^2] & \ge \delta^2 \left(1 -9\frac{C\left(\frac{L_1^2 \|K\|_2^2}{32} \vee \frac{L_2^2 \|K\|_2^4}{32}\right)\psi_n + \log{2}}{c\psi_n}\right) \\
& = \delta^2 \left(1 -9\frac{C\left(\frac{L_1^2 \|K\|_2^2}{32} \vee \frac{L_2^2 \|K\|_2^4}{32}\right) + \frac{\log{2}}{\psi_n}}{c}\right)
\end{align*}

The constants $C, c$ depends on $c_1, c_2$ via the definition of $m_1,m_2$. Choosing them appropriately and using the fact $\psi_n \uparrow \infty$ as $n \uparrow \infty$, we can make

9\frac{C\left(\frac{L_1^2 \|K\|_2^2}{32} \vee \frac{L_2^2 \|K\|_2^4}{32}\right) + \frac{\log{2}}{\psi_n}}{c} \le \frac12 \,.

which implies:

\begin{align*}
\inf_{\hat f} \sup_{f,P_{X}} \bbE_{P_{X}}[(\hat f(X) - f(X))^2] & \ge \delta^2 \\
& \ge \left(\frac{L_1^2 \|K\|_2^2}{8} \wedge \frac{L_2^2 \|K\|_2^4}{8}\right) \left( dn^{\frac{-2\beta_1}{2\beta_1 + 1}} +  \dbinom{d}{2}n^{\frac{-2\beta_2}{2\beta_2 + 2}}\right) \,.
\end{align*}

This completes the proof.

§.§ Proof of Theorem <ref> and <ref>

Recall the for any $\phi = \sum_j \phi_j + \sum_{k < l} \phi_{kl}$, we define $\|\phi\|_{n, 1} := \sum_j \|\phi_j\|_{n} + \sum_{j < k}\|\phi_{jk}\|_{n}$ and $\|\phi\|^2_n =\sum_{i=1}^{n} \phi^2(X_i)/n$. Further define $\|\phi\|_{n,lin}:=\sum_j \|\phi_j\|_{n}$ and $\|\phi\|_{n,quad}:=\sum_{k < l} \|\phi_{kl}\|_{n}$. From Theorem <ref>, there exists $\{\phi^\star_j\}_{j \in S_1}, \left\{\phi^\star_{kl}\right\}_{(j, k) \in S_2}$, $j\in S_1, (j,k) \in S_2$, such that, $\phi^\star_j \in \mathcal{F}_{NN,1}$, $\phi^\star_{jk} \in \mathcal{F}_{NN,2}$ and

\begin{align}
& \|f_{0,j} - \phi^\star_{j}\|_\infty \le C_1 (N_1L_1)^{-2\beta_1} \ \forall \ j \in S_1 \,, \label{eq:nn_approx_lin} \\
& \|f_{0,kl} - \phi^\star_{kl}\|_\infty \le C_2 (N_2L_2)^{-\beta_2} \ \forall \ (k < l) \in S_2 \,.\label{eq:nn_approx_quad}
\end{align}

and furthermore, $\phi^\star_j$ has integral $0$ and $\phi^\star_{kl}$ has $0$ marginals. For the simplicity of the rest of the proof, take $L_i$ to be of the constant order. Therefore the architecture we use here is the following: for the fitting the univariate function we use DNNs from the class $\cF(1, c_1 N_1 \log{N_1}, c_2, 1)$ and for fitting bivariate functions, we use DNNs from the class  $\cF(1, c_1 N_2 \log{N_2}, c_3, 1)$. Therefore, defining $\phi^\star = \sum_{j \in S_1} \phi^\star_j + \sum_{(j, k) \in S_2} \phi^\star_{jk}$, we have by Theorem <ref>:

\begin{align}
\bbE\left[(f_0(X)-\phi^\star(X))^2\right] & \le C_1 \left(s_1 (N_1)^{-4\beta_1}+ s_2 (N_2)^{- 2\beta_2}\right) \notag \\
& \triangleq C_1(s_1 \rho_{n, 1}^2+ s_2 \rho_{n, 2}^2)\,.\label{eq:define_nn_approx_linf}
\end{align}

As $\hat \phi$ is the minimizer of the penalized loss function among the class of neural networks, it outperforms $\phi^\star$. Setting $pen(\phi):= \lambda_{n,1} \| \phi\|_{n, lin} +\lambda_{n,2} \|\phi\|_{n, quad}$, we have:

\begin{align}
& \frac{1}{2}\|Y - \hat \phi\|_n^2 + pen(\hat\phi) \le  \frac{1}{2}\|Y - \phi^\star\|_n^2 + pen(\phi^\star) \notag \\
% \implies &  \frac{1}{2}\|\hat \phi - \phi^\star\|_n^2 + pen(\hat\phi) \le  \langle Y - \phi^\star, \hat \phi - \phi^\star\rangle_n + pen(\phi^\star) \\
\implies &  \frac{1}{2}\|\hat \phi - \phi^\star\|_n^2 + pen(\hat\phi) \le  \langle f_0 - \phi^\star, \hat \phi - \phi^\star\rangle_n + \langle \eps , \hat \phi - \phi^\star\rangle_n + pen(\phi^\star) \notag \\
\implies & \frac{1}{2}\|\hat \phi - \phi^\star\|_n^2 + pen(\hat\phi)\le \|f_0 - \phi^\star\|_n^2 + \frac14 \|\hat \phi - \phi^\star\|_n^2 + \langle \eps , \hat \phi - \phi^\star\rangle_n + pen(\phi^\star) \notag \\
\implies & \frac{1}{4}\|\hat \phi - \phi^\star\|_n^2 + pen(\hat\phi) \le \|f_0 - \phi^\star\|_n^2 + \langle \eps , \hat \phi - \phi^\star\rangle_n +pen(\phi^\star). \label{eq:fan1}
\end{align}

First we bound the empirical error $\|f_0 - \phi^\star\|_n^2$ in terms of its population counterpart $\|f_0 - \phi^\star\|_2^2$. Using Chebychev inequality, we have for $t > 0$:

\begin{equation}
\label{eq:chebychev}
\bbP\left(\left|\|f_0 - \phi^\star\|_n^2 - \|f_0 - \phi^\star\|_2^2  \right| > t\right) \le \frac{\bbE\left[(f_0 - \phi^\star)^4\right]}{nt^2}
\end{equation}

Expanding the fourth moment, we obtain:

\begin{align*}
& \bbE\left[(f_0 - \phi^\star)^4\right] \le 8\left\{\bbE\left[\left(\sum_{j \in S_1} (f_{0,j} - \phi^\star_j)\right)^4\right] + \bbE\left[\left(\sum_{(j, k) \in S_2} (f_{0,jk} - \phi^\star_{jk})\right)^4\right]\right\} \\
& \le 8p_{\max}\left\{\int_{[0, 1]^d}\left(\sum_{j \in S_1} (f_{0,j}(x_j) - \phi^\star_j(x_j))\right)^4 \ d\bx + \int_{[0, 1]^d}\left(\sum_{(j, k) \in S_2} (f_{0,jk}(x_j, x_k) - \phi^\star_{jk}(x_j, x_k))\right)^4 \ d\bx \right\} \\
& \triangleq 8p_{\max}(T_1 + T_2) \,.
% & \le p_{\max}\int_{[0, 1]^S} (f_0(x) - \phi^\star(x))^4 \ dx \\
% & = p_{\max}\sum_{j \in S_1} \int_0^1 (f_{0, j}(x_j) - \phi^\star_{n, j}(x_j))^4 \ dx_j + \sum_{(j,k) \in S_2} \iint_{[0,1]^2} (f_{0, jk}(x_j,x_k) - \phi^\star_{n, jk}(x_j,x_k))^4 \ dx_j \ dx_k
% %+ \sum_{j \neq k} \int_0^1\int_0^1 (f_{0, j}(x_j) - \phi^\star_{n, j}(x_j))^2(f_{0, k}(x_k) - \phi^\star_{n, k}(x_k))^2 \ dx_j \ dx_k
% \\
% & \le C \left(s_1 (NL)^{-8\alpha} +s_2 (NL)^{-4\alpha}\right) \,.
\end{align*}

Now we analyze $T_1$ and $T_2$ separately. For $T_1$:

\begin{align*}
T_1 & = \int_{[0, 1]^d}\left(\sum_{j \in S_1} (f_{0,j}(x_j) - \phi^\star_j(x_j))\right)^4 \ d\bx \\
& = \sum_{j \in S_1} \int_0^1 (f_{0,j}(x_j) - \phi^\star_j(x_j))^4 \ dx_j + \sum_{j \neq k \in S_1} \int_0^1 (f_{0,j}(x_j) - \phi^\star_j(x_j))^2 \ dx_j\int_0^1 (f_{0,j}(x_k) - \phi^\star_k(x_k))^2 \ dx_k \\
& \le C^4 s^2_1 (N_1L_1)^{-8\beta} \,.
\end{align*}

where the last line follows from the fact that $\|f_{0, j} - \phi^\star_j\|_\infty \le C(N_1L_1)^{-2\beta}$.
Similar analysis for $T_2$ yields:

\begin{align*}
T_2 & = \int_{[0, 1]^d}\left(\sum_{(j, k) \in S_2} (f_{0,jk}(x_j, x_k) - \phi^\star_{jk}(x_j, x_k))\right)^4 \ d\bx
% & = \sum_{(j, k) \in S_2} \int_0^1\int_0^1 (f_{0,jk}(x_j, x_k) - \phi^\star_{jk}(x_j, x_k))^4 \ dx_j \ dx_k \\
% &  \hspace{-0.3in}+ \sum_{(j, k) \neq (l, m) \in S_1} \int_{[0, 1]^4}(f_{0,jk}(x_j, x_k) - \phi^\star_{jk}(x_j, x_k))^2(f_{0,lm}(x_l, x_m) - \phi^\star_{lm}(x_l, x_m))^2  \ dx_j \ dx_k\ dx_l \ dx_m \\
% & + \hspace{-0.3in} \sum_{(j_1, k_1) \neq (j_2, k_2) \neq (j_3, k_3) \in S_1} \int_{[0, 1]^4}(f_{0,jk}(x_j, x_k) - \phi^\star_{jk}(x_j, x_k))^2(f_{0,lm}(x_l, x_m) - \phi^\star_{lm}(x_l, x_m))^2  \ dx_j \ dx_k\ dx_l \ dx_m \\
\le C^4 s^2_2 (N_2L_2)^{-4\beta} \,.
\end{align*}

using $\|f_{0, jk} - \phi^\star_{jk}\|_\infty \le C(N_2L_2)^{-\beta}$,  repeatedly and using Assumption <ref>. Using the bounds of $T_1$ of $T_2$, we have:

\begin{align*}
\bbE\left[(f_0 - \phi^\star)^4\right] &\le 8C^4p_{\max}\left(s^2_1 (N_1L_1)^{-8\beta} + s^2_2 (N_2L_2)^{-4\beta}\right)\\
&\le 8C^4p_{\max} \Big(s^2_1\rho^4_{n,1}+s^2_2\rho^4_{n,2} \Big)\,.
\end{align*}

Hence, choosing $t_0 = C^2 \sqrt{\frac{8p_{\max} (s^2_1\rho^4_{n,1}+s^2_2\rho^4_{n,2}) \log{n}}{n}}$, we have from (<ref>):

\begin{equation}
\label{eq:chebychev_2}
\bbP\left(\left|\|f_0 - \phi^\star\|_n^2 - \|f_0 - \phi^\star\|_2^2  \right| > t_0\right) \le \frac{1}{\log{n}}
\end{equation}

Define the event $\Omega_{n, 1}$ as:

\Omega_{n, 1}  = \left\{\left|\|f_0 - \phi^\star\|_n^2 - \|f_0 - \phi^\star\|_2^2  \right| \le t_0 \right\} \,.

Noticing $t_0 \ll s_1 \rho_{n, 1}^2+ s_2 \rho_{n, 2}^2$, by (<ref>), we have
on $\Omega_{n, 1}$ that

\begin{align*}
% & \frac{1}{4}\|\hat \phi - \phi^\star\|_n^2 + pen(\hat \phi) \le \|f_0 - \phi^\star\|_n^2 + \langle \eps , \hat \phi - \phi^\star\rangle_n + pen(\phi^\star)\\
% \implies &  \frac{1}{4}\|\hat \phi - \phi^\star\|_n^2 + pen(\hat \phi) \le \|f_0 - \phi^\star\|_2^2 + C_2\rho^2_n \frac{\log{n}}{\sqrt{n}}+ \langle \eps , \hat \phi - \phi^\star\rangle_n + pen(\phi^\star) \\
\eqref{eq:fan1}  \implies &  \frac{1}{4}\|\hat \phi - \phi^\star\|_n^2 + pen(\hat \phi) \le C_3(s_1 \rho_{n, 1}^2+ s_2 \rho_{n, 2}^2)  + \langle \eps , \hat \phi - \phi^\star\rangle_n + pen(\phi^\star)
\end{align*}

In the next step, we bound the empirical error, i.e. the inner product between $\eps$ and $\hat \phi - \phi^\star$. First of all, note that, by triangle inequality:

\begin{align*}
\left|\langle \eps , \hat \phi - \phi^\star\rangle_n \right| & \le
\sum_{j=1}^d \left|\frac1n \sum_i \eps_i \left(\hat \phi_j(X_{ij}) - \phi^\star_j(X_{ij})\right)\right| \\
& \qquad \qquad + \sum_{k<l} \left|\frac1n \sum_i \eps_i \left(\hat \phi_{kl}(X_{ik}, X_{il}) - \phi^\star_{kl}(X_{ik}, X_{il})\right)\right|
\end{align*}

To bound each of the summands on the RHS, we first bound each term inside the summation. Towards that direction, we use Lemma $4$ of [Fan and Gu, 2022], which says that if $\cG_n$ is a class of functions with VC dimension $V_n$, then for any fixed $g_0 \in\cG, \eps> 0, t > 0$, with probability $\ge 1 - \log{(1/\eps)}e^{-t}$:

\begin{align*}
\frac{1}{n}\left|\sum_{i=1}^n\eps_i(g(X_i) - g_0(X_i))\right| & \le C_4\left(\|g - g_0\|_n + \eps\right)\sqrt{\frac{V_{n, 1} \log{n}}{n} + \frac{t}{n}}
% & \le C_4\left(\|g - g_0\|_n + \eps\right)\sqrt{\frac{V_{n, 1} \log{p}}{n} + \frac{t}{n}}
\end{align*}

for some universal constant $c$. Note that, $\phi_j \in \cF_{NN}(1, c_1 N_1 \log{N_1}, c_2, c_3 (N_1 \log{N_1})^2, 1)$, whose VC dim $V_{n, 1} \le C_5 N_1^2\log{N_1})$ (see Lemma <ref>).
We now apply this lemma to each of the component functions with $t = 2\log{d}$ and use a union bound to conclude:

\begin{align}\label{eq:good_event_n2_1}
\sum_{j=1}^d \left|\frac1n \sum_i \eps_i \left(\hat \phi_j(X_{ij}) - \phi^\star_j(X_{ij})\right)\right| \le  C_4\sqrt{\frac{V_{n, 1} \log{p}}{n} + \frac{2\log{d}}{n}} \left(\sum_{j=1}^d\|\hat \phi_j - \phi^\star_{n, j}\|_n + d\eps_1\right)\,.
\end{align}

The above event occurs with probability $\ge 1 - \log(1/\eps_1)e^{\log{d}-2\log{d}} = 1 - \log(1/\eps_1)/d$. A similar calculation for the bivariate components yields:

\begin{align}
\label{eq:good_event_n2_2}
& \sum_{k<l} \left|\frac1n \sum_i \eps_i \left(\hat \phi_{kl}(X_{ik}, X_{il}) - \phi^\star_{kl}(X_{ik}, X_{il})\right)\right| \notag \\
& \hspace{5em}\le C_5\sqrt{\frac{V_{n, 2} \log{n}}{n} + \frac{3\log{d}}{n}} \left(\sum_{k < l}\|\hat \phi_{kl} - \phi^\star_{kl}\|_n + \frac{d(d-1)}{2}\eps_2\right) \,,
\end{align}

with probability $\ge 1 - \log{(1/\eps_2)}e^{2\log{d}-3\log{d}} = 1 - \log{(1/\eps_2)}/d$. For the rest of the analysis, define the penalty parameters $\lambda_{n, 1}, \lambda_{n, 2}$ and $\eps_1, \eps_2$ as:

\begin{align}
\label{eq:lambda_1_val} \lambda_{n, 1} & = 2C_4\sqrt{\frac{V_{n, 1} \log{p}}{n} + \frac{2\log{d}}{n}}\\
\label{eq:lambda_2_val} \lambda_{n, 2} & = 2C_5\sqrt{\frac{V_{n, 2} \log{n}}{n} + \frac{3\log{d}}{n}}  \\
\label{eq:eps_val} \eps_i & = \frac{s_i \lambda_{n, i}}{\binom{d}{i}} \text{ for } i=1,2.
\end{align}

Combining equation (<ref>) and (<ref>) we have:

\begin{align}
\label{eq:good_event_n2_3}
\left|\langle \eps, \hat \phi - \phi^\star \rangle_n\right| & \le \frac{\lambda_{n, 1}}{2} \sum_{j=1}^p\|\hat \phi_j - \phi^\star_{n, j}\|_{n, lin} +  \frac{\lambda_{n, 2}}{2} \sum_{k < l}\|\hat \phi_{kl} - \phi^\star_{kl}\|_{n, quad} \notag \\
& \hspace{15em} + \frac{s_1\lambda_{n, 1}^2}{2} + \frac{s_2\lambda_{n, 2}^2}{2} \,.
\end{align}

Define this event (<ref>) to be $\Omega_{n, 2}$ and we will later show that our choice of $(N_1, N_2)$ ensures that $\bbP(\Omega_{n, 2}) \to 1$. On the event $\Omega_{n, 1} \cap \Omega_{n, 2}$, we have:

\begin{align}
& \frac{1}{4}\|\hat \phi - \phi^\star\|_n^2 + \lambda_{n, 1} \|\hat \phi \|_{n, lin} + \lambda_{n, 2}\|\hat \phi\|_{n, quad} \notag \\
& \le C_3 (s_1\rho_{n, 1}^2 + s_2 \rho_{n, 2}^2)  + \frac{\lambda_{n, 1}}{2} \sum_{j=1}^p\|\hat \phi_j - \phi^\star_{n, j}\|_{n, lin} \notag \\
& \qquad \qquad +  \frac{\lambda_{n, 2}}{2} \sum_{k < l}\|\hat \phi_{kl} - \phi^\star_{kl}\|_{n, quad}  + \frac{s_1\lambda_{n, 1}^2}{2} + \frac{s_2\lambda_{n, 2}^2}{2} + \lambda_{n, 1} \|\hat \phi \|_{n, lin} + \lambda_{n, 2}\|\hat \phi\|_{n, quad} \notag
\end{align}

Some simple algebra yields:

\begin{align}
\label{eq:final_bound_fixed}
& \frac{1}{4}\|\hat \phi - \phi^\star\|_n^2 + \frac{\lambda_{n, 1}}{2} \|\hat \phi \|_{n, lin} + \frac{\lambda_{n, 2}}{2}\|\hat \phi\|_{n, quad} \notag \\
& \le C_3 (s_1\rho_{n, 1}^2 + s_2 \rho_{n, 2}^2)  + \frac{3\lambda_{n, 1}}{2} \sum_{j \in S_1}\|\hat \phi_j - \phi^\star_{n, j}\|_{n, lin} +  \frac{3\lambda_{n, 2}}{2} \sum_{k < l \in S_2}\|\hat \phi_{kl} - \phi^\star_{kl}\|_{n, quad} \notag \\
& \hspace{5em} + \frac{s_1\lambda_{n, 1}^2}{2} + \frac{s_2\lambda_{n, 2}^2}{2} \,.
\end{align}

If we don't assume any Restricted Strong Convexity(RSC) on the underlying function class (i.e. in our case DNNs with some pre-specified architecture), then we can use the fact $\|\hat \phi_j - \phi^\star_{n, j}\|_n \le 2B$ for $j$ and consequently we have:

\begin{align}
\label{eq:bound_fixed_no_rse}
& \frac{1}{4}\|\hat \phi - \phi^\star\|_n^2 + \frac{\lambda_{n, 1}}{2} \|\hat \phi \|_{n, lin} + \frac{\lambda_{n, 2}}{2}\|\hat \phi\|_{n, quad} \notag \\
& \le C_3 (s_1\rho_{n, 1}^2 + s_2 \rho_{n, 2}^2)  + \frac{3}{2}(s_1\lambda_{n, 1} + s_2\lambda_{n, 2}) + \frac{s_1\lambda_{n, 1}^2}{2} + \frac{s_2\lambda_{n, 2}^2}{2} \,.
\end{align}

This completes the proof of Theorem <ref>. For the second part, we assume to have RSC condition
Assumption (<ref>). From (<ref>), we obtain:

\begin{align*}
& \frac{1}{4}\|\hat \phi - \phi^\star\|_n^2 + \frac{\lambda_{n, 1}}{2} \|\hat \phi \|_{n, lin} + \frac{\lambda_{n, 2}}{2}\|\hat \phi\|_{n, quad} \notag \\
& \le C_3 (s_1\rho_{n, 1}^2 + s_2 \rho_{n, 2}^2)  + \frac{3\lambda_{n, 1}}{2} \sum_{j \in S_1}\|\hat \phi_j - \phi^\star_{n, j}\|_{n, lin} +  \frac{3\lambda_{n, 2}}{2} \sum_{k < l \in S_2}\|\hat \phi_{kl} - \phi^\star_{kl}\|_{n, quad} \notag \\
& \hspace{5em} + \frac{s_1\lambda_{n, 1}^2}{2} + \frac{s_2\lambda_{n, 2}^2}{2} \\
& \le C_3(s_1(\rho_{n, 1}^2 + \lambda_{n, 1}^2) + s_2(\rho_{n, 2}^2 + \lambda_{n, 2}^2)) + \frac{3}{2}\lambda_{n,1}\sqrt{s_1}\sqrt{\sum_{j \in S_1}\|\hat \phi_j - \phi^\star_{n, j}\|^2_{n, lin}} \\
& \qquad \qquad \qquad + \frac{3}{2} \lambda_{n, 2}\sqrt{s_2}\sqrt{\sum_{k < l \in S_2}\|\hat \phi_{kl} - \phi^\star_{kl}\|^2_{n, quad}} \\
& \le C_3 (s_1(\rho_{n, 1}^2 + \lambda_{n, 1}^2) + s_2(\rho_{n, 2}^2 + \lambda_{n, 2}^2)) + \frac{9s_1\lambda_{n, 1}^2}{2\kappa_1^2} + \frac{\kappa_1^2}{8}\sum_{j \in S_1}\|\hat \phi_j - \phi^\star_{n, j}\|^2_{n, lin} \\
& \qquad \qquad \qquad + \frac{9s_2\lambda_{n, 2}^2}{2\kappa_2^2} + \frac{\kappa_2^2}{8}\sum_{k < l \in S_2}\|\hat \phi_{kl} - \phi^\star_{kl}\|^2_{n, quad} \\
& \le C_3(s_1(\rho_{n, 1}^2 + \lambda_{n, 1}^2) + s_2(\rho_{n, 2}^2 + \lambda_{n, 2}^2)) +  \frac{9s_2\lambda_{n, 2}^2}{2\kappa_2^2} + \frac{9s_1\lambda_{n, 1}^2}{2\kappa_1^2} + \frac{1}{8}\|\hat \phi - \phi^\star\|_n^2
\end{align*}

This implies:

\begin{align*}
& \frac{1}{8}\|\hat \phi - \phi^\star\|_n^2 + \frac{\lambda_{n, 1}}{2} \sum_{j \in S_1^c}\|\hat \phi_j \|_{n} + \frac{\lambda_{n, 2}}{2}\sum_{(j < k) \in S_2^c}\|\hat \phi_{jk}\|_{n}
% & \qquad \qquad \le 2(s_1(\rho_{n, 1}^2 + \lambda_{n, 1}^2) + s_2(\rho_{n, 2}^2 + \lambda_{n, 2}^2)) +  \frac{9s_2\lambda_{n, 2}^2}{2\kappa_2^2} + \frac{9s_1\lambda_{n, 1}^2}{2\kappa_1^2} \\
& \le C_6\left(s_1(\rho_{n, 1}^2 + \lambda_{n, 1}^2) + s_2(\rho_{n, 2}^2 + \lambda_{n, 2}^2)\right)
\end{align*}

Now we select the optimal value of $N_i$'s by balancing $\rho_{n, i}$ and $\lambda_{n, i}$. First take $i = 1$. Recall that we have:

\rho_{n, 1}^2 = N_1^{-4\beta_1}, \text{   and   } \lambda^2_{n, 1} = 4C_4^2 \left(\frac{V_{n, 1}\log{n}}{n} + \frac{2\log{d}}{n}\right) \le C_7 \left(\frac{N_1^2\log^3{N_1}\log{n}}{n}  + \frac{\log{d}}{n}\right)

Choosing $N_1 = \lfloor n^{\frac{1}{2(2\beta_1 + 1)}}\rfloor$ yields:

\begin{align}
\label{eq:rho_val_final}
\rho^2_{n, 1} & \sim n^{-\frac{2\beta_1}{2\beta_1 + 1}} \,, \\
\label{eq:lambda_val_final}
\lambda^2_{n, 1} & \sim n^{-\frac{2\beta_1}{2\beta_1 + 1}}\log^4{n} + \frac{\log{d}}{n} \,.
\end{align}

In particular we have:

\rho_{n, 1}^2 + \lambda^2_{n, 1} \le C_8 \left(n^{-\frac{2\beta_1}{2\beta_1 + 1}}\log^4{n} + \frac{\log{d}}{n}\right) \,.

Similar calculation for $i = 2$ implies that choosing $N_2 = \lfloor n^{\frac{1}{2(\beta_2 + 1)}}\rfloor$ we have:

\rho_{n, 2}^2 + \lambda^2_{n, 2} \le C_9 \left(n^{-\frac{\beta_2}{\beta_2 + 1}}\log^4{n} + \frac{\log{d}}{n}\right) \,.

Therefore, the above choice of penalty parameters yields:

\begin{align*}
& \frac{1}{8}\|\hat \phi - \phi^\star\|_n^2 + \frac{\lambda_{n, 1}}{2} \sum_{j \in S_1^c}\|\hat \phi_j \|_{n} + \frac{\lambda_{n, 2}}{2}\sum_{(j < k) \in S_2^c}\|\hat \phi_{jk}\|_{n} \\
& \qquad \qquad \le C_{10}\left(s_1\left(n^{-\frac{2\beta_1}{2\beta_1 + 1}}\log^4{n} + \frac{\log{d}}{n}\right) + s_2\left(n^{-\frac{\beta_2}{\beta_2 + 1}}\log^4{n} + \frac{\log{d}}{n}\right)\right) \,.
\end{align*}

This completes the proof.

§.§ Proof of Theorem <ref>