Independence-Encouraging Subsampling for Nonparametric Additive Models

TL;DR

This paper introduces independence-encouraging subsampling (IES), a novel method for selecting representative subsamples from large datasets to efficiently fit nonparametric additive models while maintaining statistical properties.

Contribution

The paper proposes IES, a subsampling technique inspired by orthogonal arrays, to improve the efficiency of additive model fitting on large datasets without sacrificing theoretical guarantees.

Findings

IES subsample converges to an orthogonal array

Backfitting over IES subsample converges to a unique solution

Numerical experiments validate the effectiveness of IES

Abstract

The additive model is a popular nonparametric regression method due to its ability to retain modeling flexibility while avoiding the curse of dimensionality. The backfitting algorithm is an intuitive and widely used numerical approach for fitting additive models. However, its application to large datasets may incur a high computational cost and is thus infeasible in practice. To address this problem, we propose a novel approach called independence-encouraging subsampling (IES) to select a subsample from big data for training additive models. Inspired by the minimax optimality of an orthogonal array (OA) due to its pairwise independent predictors and uniform coverage for the range of each predictor, the IES approach selects a subsample that approximates an OA to achieve the minimax optimality. Our asymptotic analyses demonstrate that an IES subsample converges to an OA and that the backfitting algorithm over the subsample converges to a unique solution even if the predictors are highly dependent in the original big data. The proposed IES method is also shown to be numerically appealing via simulations and a real data application.

Tables

Table 1. Table 1 : Estimation and prediction performances of Rand and IES in the diamond price prediction data.

	Rand	IES
$\mathrm{ASE}$	$0.13$	$0.01$
$\mathrm{MEE}$	$1.51$	$0.48$
AvePredError	$0.06$	$0.06$
MaxPredError	$1.67$	$1.29$

Table 2. Table 2 : Average computation times (in seconds) spent on subsampling, CV, and model fitting. Standard deviations (SD) are in parentheses.

	Full	Rand	IES
Subsampling	$0$ $\left(0\right)$	$0.0003$ $\left(0.0000\right)$	$5.82$ $\left(0.34\right)$
CV	$8092.53$ $\left(121.99\right)$	$926.06$ $\left(20.12\right)$	$1140.04$ $\left(27.52\right)$
Fitting	$0.21$ $\left(0.11\right)$	$0.03$ $\left(0.02\right)$	$0.03$ $\left(0.01\right)$
Total	$8092.74$ $\left(121.96\right)$	$926.09$ $\left(20.12\right)$	$1145.89$ $\left(27.48\right)$

Equations

$$Y_i=m(X_i)+{ϵ}_i,i=1,\dots ,N,$$

$$\tilde{m}(x)=\sum _{i=1}^N{w}_i(x;X_1,\dots ,X_N)Y_i,$$

$$\mathcal{C}^{*}=\left\{m(x)\in C^{(2)}[0,1]\quad\big{|}\quad\max_{x}|m^{(2)}(x)|^{2}\leq\eta\right\},$$

$$R_0(x)=\frac{3}{4}15^{-1/5}{\left\{\frac{{\eta }^{1/4}{\sigma }^2}{Nf(x)}\right\}}^{4/5}\{1+o_P(1)\},$$

$$R(f)=\underset{\tilde{m}(x)\text{mbox}linear}{\min }\underset{m\in {\mathcal{C}}^{\ast },x\in [0,1]}{\sup }\mathrm{E}[(\tilde{m}(x)-m(x){)}^2\mid X_1,\dots ,X_N],$$

$$R(f)=\frac{3}{4}15^{-1/5}{\eta }^{1/5}{\left(\frac{{\sigma }^2}{Nf(x_0)}\right)}^{4/5}(1+o_p(1)),$$

$$Y_i=m({\mathbf{X}}_i)+{ϵ}_i=\mu +m_1(X_{i1})+m_2(X_{i2})+\cdots +m_p(X_{ip})+{ϵ}_i,i=1,2,\dots ,N,$$

$${\int }_0^1{m}_j(x)dx=0,j=1,\dots ,p.$$

$$m_j(X_{ij})-\mathrm{E}[m_j(X_{ij})]=\mathrm{E}[Y_i\mid X_{ij}]-\mathrm{E}[Y_i],\text{for each j=1,…,p,}$$

$$\left(\begin{array}{ccc}0& 0& 0\\ 0& 1& 1\\ 1& 0& 1\\ 1& 1& 0\end{array}\right).$$

$$X_{ij}=\frac{a_{ij}+U_{ij}}{q},\text{mbox}fori=1,\dots ,N,\text{mbox}andj=1,\dots ,p,$$

$$z_{ij}=\lfloor x_{ij}q\rfloor ,$$

$$\delta ({\mathbf{x}}_i,{\mathbf{x}}_{i^{\prime }})=\sum _{j=1}^p\text{mathbbm}1(\lfloor x_{ij}q\rfloor =\lfloor x_{i^{\prime }j}q\rfloor )=\sum _{j=1}^p\text{mathbbm}1(z_{ij}=z_{i^{\prime }j}),$$

$$L({\mathcal{X}}^{\ast })=\sum _{1\le i<i^{\prime }\le n}[\delta ({\mathbf{x}}_i^{\ast },{\mathbf{x}}_{i^{\prime }}^{\ast }){]}^2.$$

$$L({\mathcal{X}}^{\ast })\ge \frac{n}{2q^2}[np(p+q-1)-(pq{)}^2],$$

$${\mathcal{X}}_{opt}^{\ast }=\arg \underset{{\mathcal{X}}^{\ast }\subset \mathcal{X}}{\min }L({\mathcal{X}}^{\ast }).$$

$$F_n(x_1,x_2)=\frac{1}{n}\sum _{i=1}^n\text{mathbbm}1(X_{ij}^{\ast }\le x_1,X_{i{j}^{\prime }}^{\ast }\le x_2).$$

$$\underset{x_1,x_2\in [0,1]}{\sup }\left|F_n(x_1,x_2)-x_1{x}_2\right|=O_p\left(N^{-\nu /2}\right).$$

$${\mathbf{x}}_{k+1}^{\ast }$$

$$=\arg \underset{\mathbf{x}\in \mathcal{X}/{\mathcal{X}}_{(k)}^{\ast }}{\min }\left\{\sum _{1\le i<i^{\prime }\le k}[\delta ({\mathbf{x}}_i^{\ast },{\mathbf{x}}_{i^{\prime }}^{\ast }){]}^2+\sum _{1\le i\le k}[\delta ({\mathbf{x}}_i^{\ast },\mathbf{x}){]}^2\right\}$$

$$=\arg \underset{\mathbf{x}\in \mathcal{X}/{\mathcal{X}}_{(k)}^{\ast }}{\min }l(\mathbf{x}\mid {\mathcal{X}}_{(k)}^{\ast }),$$

$$l(\mathbf{x}\mid {\mathcal{X}}_{(k)}^{\ast })=\sum _{1\le i\le k}[\delta ({\mathbf{x}}_i^{\ast },\mathbf{x}){]}^2$$

$$l\left(\mathbf{x}\mid {\mathcal{X}}_{(k+1)}^{\ast }\right)=l\left(\mathbf{x}\mid {\mathcal{X}}_{(k)}^{\ast }\right)+\delta (\mathbf{x},{\mathbf{x}}_{k+1}^{\ast }{)}^2,$$

$$m(\mathbf{X})=1+\frac{8}{4+X_1}+\frac{\exp \left\{3-X_2^2\right\}}{4}+1.5\sin \left(\frac{\pi }{2}{X}_3\right),$$

$$\text{price}\approx \mu +m_1(\text{carat})+m_2(\text{depth})+m_3(\text{table})$$

$$V_{nt}(x)=\frac{1}{n}\sum _{i=1}^n\frac{1}{h_1}K\left(\frac{x_{i1}^{\ast }-x}{h_1}\right)(x_{i1}^{\ast }-x{)}^t,\text{and}{W}_{nt}(x)=\frac{1}{n}\sum _{i=1}^n\frac{1}{h_2}K\left(\frac{x_{i2}^{\ast }-x}{h_2}\right)(x_{i2}^{\ast }-x{)}^t.$$

$$[{\mathcal{S}}_1{]}_{ij}$$

$$\text{and}[{\mathcal{S}}_2{]}_{ij}$$

$$\left(\begin{array}{cc}\mathcal{I}& {\mathcal{S}}_1^{\ast }\\ {\mathcal{S}}_2^{\ast }& \mathcal{I}\end{array}\right)\left(\begin{array}{c}{\hat{\mathbf{m}}}_1\\ {\hat{\mathbf{m}}}_2\end{array}\right)=\left(\begin{array}{c}{\mathcal{S}}_1^{\ast }\\ {\mathcal{S}}_2^{\ast }\end{array}\right)\mathbf{Y},$$

$$\left(\begin{array}{c}{\hat{\mathbf{m}}}_1\\ {\hat{\mathbf{m}}}_2\end{array}\right)=\left(\begin{array}{c}\left[\mathcal{I}-{\left(\mathcal{I}-{\mathcal{S}}_1^{\ast }{\mathcal{S}}_2^{\ast }\right)}^{-1}\left(\mathcal{I}-{\mathcal{S}}_1^{\ast }\right)\right]\mathbf{Y}\\ \left[\mathcal{I}-{\left(\mathcal{I}-{\mathcal{S}}_2^{\ast }{\mathcal{S}}_1^{\ast }\right)}^{-1}\left(\mathcal{I}-{\mathcal{S}}_2^{\ast }\right)\right]\mathbf{Y}\end{array}\right).$$

Full text

Keywords: Empirical independence; Local polynomial regression; Minimax risk; Optimal design; Orthogonal array.

1 Introduction

Big data of huge sample sizes are prevalent in many disciplines such as science, engineering, and medicine. Such data may reveal important domain knowledge, but meanwhile they pose challenges to data storage and analysis. To address those challenges, subsampling has recently received increasing attention and has been intensively studied.

An optimal subsampling approach typically specifies a downstream model and carefully selects an informative subsample so that the model training on the subsample is more accurate than that on other possible subsamples. Different subsampling approaches have been developed for various parametric models. For linear regression, Ma and Sun (2015) proposed subsampling probabilities defined via leverage scores. Wang et al. (2019) investigated an information based optimal subsampling algorithm motivated by $D$-optimal experimental design. Wang et al. (2021) developed an orthogonal subsampling (OSS) method inspired by the universal optimality of orthogonal array (OA) for linear regression. Subsampling methods for other parametric models are also extensively studied, such as Wang et al. (2018) and Han et al. (2020) for logistic regressions, Wang and Ma (2021) for quantile regression, and Ai et al. (2021) for generalized linear models. Despite their optimality in some sense for fitting specific parametric models, the usage of those methods can be hindered by strong model assumptions that may not hold in big data problems. See Fan et al. (2014) for a detailed discussion. To this end, Meng et al. (2021) proposed an algorithm, called LowCon, to select a space-filling subsample which is shown to be robust when a linear model is misspecified. Researchers have also looked into nonparametric settings with less stringent model assumptions. For example, Meng et al. (2020) showed the superiority of a space-filling subsample for multivariate smoothing splines; Yang et al. (2017) applied tensor sketching to accelerate kernel ridge regression; Zhao et al. (2018) and He and Hung (2022) considered design-based subsampling for Gaussian process modeling; Shi and Tang (2021) considered model-robust subdata selection. Other methods include continuous distribution compression (Mak and Joseph, 2018) and supervised data compression (Joseph and Mak, 2021).

The nonparametric additive model (Hastie and Tibshirani, 1986) has been widely used in practice because of its interpretability and flexibility (e.g., Walker and Wright, 2002; Hwang et al., 2009; Liutkus et al., 2014). It avoids the “curse of dimensionality” which impedes the implementation of fully nonparametric models with multiple predictors. However, fitting an additive model may still be computationally expensive when the sample size is huge. For example, if the backfitting algorithm (Breiman and Friedman, 1985; Buja et al., 1989) combined with local polynomial smoothing is used to fit an additive model on a data set with $N$ observations of $p$ predictors, where $p\ll N$, the time complexity is $O(N^{2})$ per backfitting iteration. If the bandwidth is selected via cross-validation, then the complexity would become $O(N^{2})$ per bandwidth grid evaluation. Therefore, the practicality of additive models is hindered for large data.

We propose an independence-encouraging subsampling (IES) method for fitting an additive model with big data. Akin to the OSS (Wang et al., 2021), the IES is inspired by the robustness and optimality of OA for experimental design and data collection (Cheng, 1980; Taguchi and Clausing, 1990). Nevertheless, existing results for OAs focus on their optimality for identifying main effects and interactions via linear regression. We first derive theoretical results on the minimax optimality of random OAs for nonparametric additive models and then develop the IES method to select a subsample that approximates a random OA. The merits of IES are three-fold. Firstly, it is fast and easy to implement. The computation of selecting a subsample and training a nonparametric additive model on the subsample is significantly faster than training the model on the large full data. Secondly, our theoretical analyses show that an IES subsample converges to a random OA whose predictors achieve marginal uniformity and pairwise independence. This substantially benefits the backfitting algorithm, the most popular numerical approach to fit additive models. A well-known sufficient condition for local polynomial backfitting estimator to converge is the “near independence” between predictors (Opsomer and Ruppert, 1997). Since the predictors are empirically independent in the selected subsample, the nbackfitting algorithm will converge to a unique solution even if the predictors are highly dependent in the original big data. Lastly, the IES approach is numerically shown to be superior to existing subsampling methods for fitting additive models and robust against certain model misspecifications.

The remainder of this paper proceeds as follows. Section 2 derives the minimax optimal sampling plan for additive models. Section 3 introduces random OAs and their properties. Section 4 proposes the IES subsampling approach, and develops some asymptotic theories. Section 5 provides a fast implementation algorithm for IES. Sections 6 and 7 present simulations and a real data example, respectively. Discussion in Section 8 concludes this paper. Technical proofs are provided in the Supplementary Materials. R code is publicly available at https://github.com/.....

2 Minimax Optimal Sampling Plan

In this section, we introduce the minimax optimal sampling plan for univariate nonparametric regression and then extend it to additive models.

2.1

Optimal sampling for univariate nonparametric regression

We first consider univariate nonparametric regression for independent and identically distributed (i.i.d.) data:

[TABLE]

where for the $i$-th subject, $i=1,\ldots,N$, $X_{i}$ is the univariate continuous predictor, $Y_{i}$ is the response, $\epsilon_{i}$ is the random error, and $m(x)=E(Y_{i}\mid X_{i}=x)$ is the regression function. The support of the predictor is assumed compact and hereafter $[0,1]$ without loss of generality. It is also assumed that $\epsilon_{i}$ are independent of the predictors, $\mathrm{E}(\epsilon_{i})=0$, and $\mathrm{Var}(\epsilon_{i})=\sigma^{2}$.

The literature of univariate nonparametric regression (e.g., Chapter 5 of Wasserman, 2006) favors linear smoothers of the form

[TABLE]

where $w_{i}$ is a data-dependent weight function. Among them, the local linear estimator is a popular option. Fan (1992) showed that under mild conditions, the local linear estimator for (1) asymptotically achieves a minimax risk on the mean squared error (MSE), where the minimum is taken over all linear smoothers and the maximum is taken over all $m(\cdot)$ in

[TABLE]

with $C^{2}[0,1]$ denoting the set of functions whose second derivatives are continuous. For any $x\in[0,1]$, the minimax risk is

[TABLE]

where $f$ is the density of the predictor distribution.

The $R_{0}(x)$ may still be large for the region with a small $f(x)$. We hope that an estimator is “robust” for all $m(\cdot)$ in $\mathcal{C}^{*}$ and all $x\in[0,1]$, in the sense that the estimator performs well even in the worst scenario. Therefore, we seek a sampling regime, or equivalently a design density $f$, that minimizes the following minimax risk:

[TABLE]

where $\min_{\tilde{m}(x)\mbox{\footnotesize~{}linear}}$ takes the minimum over all linear smoothers in (2), and $\mathcal{C}^{*}$ is defined in (3). The following result calculates the $R(f)$ in (4) and provides the optimal $f$ that minimizes $R(f)$. Denote $[a]_{+}=\max\{0,a\}.$

Theorem 1**.**

Suppose that $f(x)$ is bounded away from zero and infinity. Let $f(x_{0})=\min_{x\in[0,1]}f(x)$. The minimax risk in (4) is given by

[TABLE]

which is achieved by the local linear regression estimator with the Epanechnikov kernel $K_{0}(u)=3[1-u^{2}]_{+}/4$ and bandwidth $h_{0}=\left\{15\sigma^{2}/[N\eta f(x_{0})]\right\}^{1/5}.$ The optimal design density $f$ that minimizes $R(f)$ in (5) is the uniform density, that is, $f(x)=1$ for all $x\in[0,1]$.

2.2

Optimal sampling for additive models

We now consider an additive model for i.i.d. data:

[TABLE]

where $\mathbf{X}_{i}=(X_{i1},\dots,X_{ip})$ contains $p$ predictors, $Y_{i}$ is the response, $m(x)=E(Y_{i}|\mathbf{X}_{i}=x)$ is the regression function, $\mu$ is a constant, $m_{j}(x)$ is the component function for the $j$-th predictor assumed to be smooth, and $\epsilon_{i}$’s are random errors. Again, the support of each predictor is assumed $[0,1]$ without loss of generality, and $\epsilon_{i}$ is independent of the predictors with $\mathrm{E}(\epsilon_{i})=0$ and $\mathrm{Var}(\epsilon_{i})=\sigma^{2}$. Moreover, the following condition is imposed for identifiability:

[TABLE]

The backfitting algorithm (e.g., Breiman and Friedman, 1985; Buja et al., 1989) is a popular, intuitive, and easy-to-implement numerical approach for fitting additive models. The algorithm updates each component function estimator alternately and iteratively. At each iteration, a one-dimensional smoother, e.g., the local linear smoother, is applied to regress the residual on one predictor to update its corresponding component function estimate, where the residual is obtained by subtracting all the other component functions’ estimates from the response. The asymptotic properties of the backfitting algorithm have been studied by Opsomer and Ruppert (1997) and Opsomer (2000). Their results also indicate that the convergence of the backfitting algorithm is not theoretically guaranteed if some predictors are highly dependent.

On the contrary, if all predictors are pairwise independent, (6) implies that

[TABLE]

where the left-hand side is a centered component function and the right-hand side suggests a univariate regression of the response on the $j$-th predictor. Hence, pairwise independence separates the additive modeling problem to $p$ one-dimensional estimations, so no iteration is required. In fact, as shown in Opsomer and Ruppert (1997), “near independence” between predictors can ensure the local-polynomial-based backfitting algorithm to converge. Therefore, inspired by Theorem 1, we recommend sampling predictors independently and uniformly to achieve the minimax optimality for each component function estimation. By Theorem 3.1 of Opsomer (2000), marginal uniformity is also optimal in minimizing the conditional variance of each local polynomial-based backfitted component function estimator over all possible designs.

When selecting a subsample from large data, since the data may have highly dependent predictors and follow an arbitrary distribution, obtaining a subsample with independently and uniformly distributed predictors (at the population level) is typically impossible. However, we can seek empirical independence and uniformity for predictors in the subsample, and this can be achieved via random OA.

3 Introduction to OA

An OA of strength $t$, denoted by OA$(N,p,q,t)$, is an $N\times p$ matrix with entries of $q$ levels indexed by $\{0,1,2,\dots,q-1\}$, arranged in such a way that all level combinations occur equally often in any $t$ columns (Hedayat et al., 1999). Such equal frequency of level combinations is called combinatorial orthogonality. The following matrix, as an example, is an OA$(4,3,2,2)$, any two columns of which consist of $(0,0)$, $(0,1)$, $(1,0)$, and $(1,1)$ exactly once:

[TABLE]

In this paper, OAs mentioned are assumed to have strength 2 unless otherwise specified.

OAs have been extensively used as fractional factorial designs because they allow uncorrelated estimation of main effects through linear regression (Wu and Hamada, 2011; Mukerjee and Wu, 2006; Wang and Xu, 2022). Cheng (1980) showed that an OA on $q$ levels is universally optimal, i.e., optimal under a wide variety of criteria that include $D$- and $A$-optimality, among all $q$-level factorial designs for studying main effects.

We now extend the superiority of OAs for establishing nonparametric additive models. Consider the sampling distribution of the column variables $A_{j}$ in an OA. We have $P(A_{j}=a)=q^{-1}$ and $P(A_{j}=a,A_{j^{\prime}}=a^{\prime})=P(A_{j}=a)P(A_{j^{\prime}}=a^{\prime})=q^{-2}$ for all $a,a^{\prime}\in\{0,1,2,\dots,q-1\}.$ Therefore, any column variable in an OA follows a discrete uniform distribution, and any pair of column variables are independent. We next provide a sampling scheme to draw data from $[0,1]^{p}$ that carry over the uniformity and variable independence of an OA.

Definition 1**.**

Given an OA$(N,p,q,2)$, denoted by $\mathcal{A}=(a_{ij})$ for $i=1,\ldots,N$ and $j=1,\ldots,p$, a random OA $(X_{ij})$ is given by

[TABLE]

where the $U_{ij}$’s are independent uniform random variables on $[0,1]$.

A random OA can be understood as a two-step sampling procedure. Firstly, partition the cube $[0,1]^{p}$ into $q^{p}$ equal-sized cells (subcubes with each side of length $q^{-1}$) and select the $n$ cells specified by the rows of $\mathcal{A}$. The $i$th row of $\mathcal{A}$ specifies the cell $\Pi_{j=1}^{p}[a_{ij}/q,(a_{ij}+1)/q)$. Secondly, randomly draw a point from each selected cell. Figure 1 illustrates the four selected cells according to (8). For any two columns, the projection of selected cells covers the whole face. Therefore, the randomly sampled points from those cells uniformly cover any two-dimensional subspace. Such a sampling scheme was also studied in Owen (1992) to obtain a better approximation of integration than Monte Carlo sampling.

Lemma 1**.**

For a random OA, the cumulative distribution on each column is given by $F(x_{1})=x_{1}$, and on any pair of columns is given by $F(x_{1},x_{2})=x_{1}x_{2}$.

Lemma 1 claims both uniformity and pairwise independence between column variables in a random OA, which are inherited from its combinatorial orthogonality and are the exact properties we seek for the optimal training data for additive models. It should be noted that for an OA$(N,p,q,2)$ to exist, the number of rows has to be a multiple of $q^{2}$, that is, $N=\lambda q^{2}$ for some positive integer $\lambda$. Abundant methods have been proposed to generate OAs, and we relegate a summary of their wide availability and generating methods to Appendix A.

4 Independence-Encouraging Subsampling (IES)

Let $(\mathbf{x}_{1},y_{1}),\dots,(\mathbf{x}_{N},y_{N})$ denote the full data with $N$ observations, where $\mathbf{x}_{i}=(x_{i1},\ldots,x_{ip})$ are observations of $p$ predictors and $y_{i}$ is the corresponding response. We consider taking a subsample of size $n$, denoted as $(\mathbf{x}_{1}^{*},y_{1}^{*}),\dots,(\mathbf{x}_{n}^{*},y_{n}^{*})$. Based on the previous discussion, our goal is to encourage empirical uniformity and pairwise independence of predictors in the subsample, and this can be achieved by finding a subsample whose design matrix approximates a random OA.

An intuitive approach is to choose an existing OA with $n$ rows and randomly select a data point in each cell specified by the OA. This approach has two possible limitations. First, for an OA$(n,p,q,2)$ to exist, the number of rows has to be a multiple of $q^{2}$, meaning that this approach is possible only when $n=\lambda q^{2}$ for some positive integer $\lambda$. Second, even if $n$ is a multiple of $q^{2}$, the full data may not fit an arbitrarily chosen OA, that is, many cells of the OA may be empty and do not contain any data points.

The proposed IES method selects a subsample by directly minimizing a discrepancy function that measures its deviation from an OA. As a result, the subsample size is not restricted to be a multiple of a square number, and the selected subsample approximates an OA that is the best compatible with the data.

4.1 The IES approach

For a full data with design matrix $\mathcal{X}=(x_{ij})$ and a prespecified integer $q$, define the membership matrix as $\mathcal{Z}=(z_{ij})$, where

[TABLE]

for $i=1,2,\dots,N$, and $j=1,2,\dots,p.$ Clearly $z_{ij}\in\{0,1,2,\dots,q-1\}$. Our goal is to search for a subsample whose design matrix $\mathcal{X}^{*}$ has an OA membership matrix. For any two observations with $\mathbf{x}_{i}$ and $\mathbf{x}_{i^{\prime}}$, define

[TABLE]

where $\mathbbm{1}{(z_{ij}=z_{i^{\prime}j})}$ is the indicator function that equals 1 if $z_{ij}=z_{i^{\prime}j}$ and 0 otherwise. Here, $\delta(\mathbf{x}_{i},\mathbf{x}_{i^{\prime}})$ counts the membership coincidence between elements of $\mathbf{z}_{i}$ and $\mathbf{z}_{i^{\prime}}$, and thus measures the similarity between $\mathbf{x}_{i}$ and $\mathbf{x}_{i^{\prime}}$. For a subsample with design matrix $\mathcal{X}^{*}=(\mathbf{x}^{*}_{1},\dots,\mathbf{x}^{*}_{n})^{T}$, define

[TABLE]

Clearly, $L(\mathcal{X}^{*})$ measures the overall similarity between all data points in $\mathcal{X}^{*}$. The following theorem shows that $L(\mathcal{X}^{*})$ also measures the discrepancy between $\mathcal{X}^{*}$ and an OA.

Theorem 2**.**

For any $\mathcal{X}^{*}$ with $n$ rows,

[TABLE]

and the lower bound is achieved if and only if $\mathcal{Z}$, the membership matrix of $\mathcal{X}^{*}$, is an OA$(n,p,q,2)$.

Theorem 2 shows that $L(\mathcal{X}^{*})$ has a lower bound which is attained if and only if the membership matrix of $\mathcal{X}^{*}$ forms an OA. In this sense, $L(\mathcal{X}^{*})$ can be viewed as a metric on the discrepancy between $\mathcal{X}^{*}$ and a realization of a random OA. Therefore, we propose the IES method, which solves the optimization problem:

[TABLE]

The IES subsample is $\{\mathcal{X}^{*}_{opt},\mathbf{y}^{*}_{opt}\}$, where $\mathbf{y}^{*}_{opt}$ is the corresponding response vector.

The optimization in (10) does not impose any restriction on $n$. When $n=\lambda q^{2}$ and an OA$(n,p,q,2)$ exists, we obtain a subsample from (10) with an OA membership matrix. Otherwise, we obtain a subsample that approximates the combinatorial orthogonality in an OA. We can extend Theorem 2 to a more general setting of $n$ for which an OA$(n,p,q,2)$ may not exist, which confirms that the optimization in (10) best approximates an OA for a general setting of $n$. The presentation of the result requires tedious notations and concepts, so we relegate the details to Lemma SLABEL:subsampwoa in Supplementary Material.

We next investigate the asymptotic properties of an IES subsample selected by (10), under the following assumptions.

Assumption 1**.**

The probability density function that generates the design matrix of the full data is compactly supported and bounded away from zero and infinity.

Assumption 2**.**

There exists some fixed positive integer $\lambda$ such that $n-\lambda q^{2}=O(q)$, and an OA$(q^{2},p+1,q,2)$ exists.

Assumption 3**.**

The subsample size $n$ goes to $\infty$ at the rate of $O(N^{\nu})$ for some $\nu\in\left(0,2/p\right)$.

Assumption 1 ensures that the full data asymptotically cover the design region as the size $N$ increases. Assumption 2 indicates again that the IES does not require $n=\lambda q^{2}$. The requirement of the existence of OA$(q^{2},p+1,q,2)$ is weak, as discussed in Appendix A, especially considering that we can set $q$ to be much bigger than $p$. Assumption 3 requires that $n$ does not grow faster than $N^{2/p}$, which is commonly the case in the setting of big subsampling.

Theorem 3**.**

Define the induced joint cumulative distribution function on any two columns of $\mathcal{X}^{*}_{opt}$, $X_{j}^{*}$ and $X_{j^{\prime}}^{*}$, as

[TABLE]

Then under Assumptions 1-3, we have

[TABLE]

Theorem 3 shows that asymptotically the solution to (10) achieves pairwise independence and uniformity, leading to a desired subsample for additive models. The convergence rate depends on $\nu$ in Assumption 3. A bigger $\nu$ indicates a larger subsample size and results in a faster convergence to the uniform distribution. We can relax Assumption 2 to a more general setting of $n$ with $n-\lambda q^{2}=O(q^{\gamma})$ for some $\gamma\in(0,2)$. The case of $\gamma\leq 1$ is equivalent to Assumption 2, and for $\gamma>1$, $F_{n}(x_{1},x_{2})$ still converges to uniformity but at a slower rate; see the proof of Theorem 3 in the Supplementary Materials for details.

4.2 Additive Modeling on IES Subsamples

After obtaining the subsample $\{\mathcal{X}^{*}_{opt},\mathbf{y}^{*}_{opt}\}$ from (10), we fit an additive model on this subsample. Since the predictors in the subsample cannot be guaranteed to be perfectly independent, we propose to estimate each component function via the backfitting algorithm (Breiman and Friedman, 1985). Motivated by Theorem 1, we apply local linear smoothers in each backfitting step.

When there are two predictors, i.e., $p=2$, we can prove the convergence of the backfitting algorithm on the subsample $\{\mathcal{X}^{*}_{opt},\mathbf{y}^{*}_{opt}\}$. We need the following assumptions in addition to Assumptions 1–3.

Assumption 4**.**

The kernel function $K$ is a symmetric density function compactly supported on $[-1,1]$. Moreover, $K$ is $M$-Lipschitz for some constant $M>0$, i.e., $|K(u)-K(v)|\leq M|u-v|$ for any $u,v\in[-1,1]$.

Assumption 5**.**

As the size of the subsample $n\rightarrow\infty$, the bandwidth $h_{j}\rightarrow 0$ and $nh_{j}^{4}\rightarrow\infty$ for $j=1,2$.

Both Assumptions 4 and 5 will be used in the proof of Theorem 4 to control certain numerical integration errors. Assumption 4 on the kernel function is commonly adopted by kernel-smoothing-based additive modeling methods (e.g., Opsomer and Ruppert, 1997; Zhang et al., 2013) and can be satisfied by popular kernels, e.g., the Epanechnikov kernel. Assumption 5 on bandwidths is mild and can be satisfied if each $h_{j}$ takes the optimal order $n^{-1/5}$ as in the literature of local polynomial smoothing (e.g., Fan and Gijbels, 1996).

Theorem 4**.**

*Under Assumptions 1-5, when $p=2$, the backfitting algorithm on the subsample $\{\mathcal{X}^{*}_{opt},\mathbf{y}^{*}_{opt}\}$ converges to a unique solution with probability approaching one as $N\rightarrow\infty$. *

The expression of the unique solution involves more tedious notations and can be unwieldy in practice. To save space, we defer the details to Appendix B. Substantially different from the result by Opsomer and Ruppert (1997), Theorem 4 does not require a weak dependency between the two predictors in the population; even if the population dependency between the predictors is high, they are almost independent in $\mathcal{X}^{*}_{opt}$ as guaranteed by Theorem 3, so the backfitting procedure on the subsample can converge asymptotically. Another critical distinction between Theorem 4 and Opsomer and Ruppert (1997) is that the latter handles independent observations while observations in $\mathcal{X}^{*}_{opt}$ are dependent.

When $p\geq 3$, theoretical convergence for the backfitting procedure on an IES subsample is unknown and will be deferred for future work. Nevertheless, it always converges numerically in our simulation studies and real data application in Sections 6 and 7.

5 Practical Implementation of IES

The optimization problem in (10) is computationally expensive to solve. An exhausted search requires evaluating the quantity $L(\mathcal{X}^{*})$ on $\binom{N}{n}$ possible subsamples, which is prohibitive for even a moderate data size. To improve the efficiency, we propose a sequential IES implementation which selects subsample points sequentially. We start with a randomly selected point $(\mathbf{x}_{1}^{*},y_{1}^{*})$. Denote the subsample design matrix with $k$ points as $\mathcal{X}^{*}_{(k)}=(\mathbf{x}^{*}_{1},\dots,\mathbf{x}^{*}_{k})^{T}$ for $k\in\{1,\ldots,n-1\}$. The $(k+1)$th subsample point is then selected as

[TABLE]

where

[TABLE]

measures the similarity between $\mathbf{x}$ and $\mathcal{X}^{*}_{(k)}$, and the selected $\mathbf{x}^{*}_{k+1}$ is the least similar point to $\mathcal{X}^{*}_{(k)}$. If there are multiple minimizers, $\mathbf{x}^{*}_{k+1}$ is randomly selected among them. After choosing $\mathbf{x}^{*}_{k+1}$, we update $l(\cdot)$ for $\mathbf{x}\in\mathcal{X}/\mathcal{X}^{*}_{(k+1)}$ via

[TABLE]

so the computational complexity of selecting one point is $O(Np)$.

Algorithm 1 outlines the detailed steps of the sequential IES implementation. In our numerical results in Sections 6 and 7, the backfitting algorithm uses the local linear smoothing and is conducted via the R package gam (Hastie, 2015). The hyperparameter $q$ can be any not-to-small integer, and we find that an integer greater than 10 would be adequate. Also, setting $q$ at a prime power may provide more stable numerical performance because of the better OA approximation and combinatorial orthogonality (details in Appendix A). Therefore, we recommend choosing a prime power $q$ which is close to $\sqrt{n/\lambda}$ for some positive integer $\lambda$. In our simulation and real data studies where $n=1000$ and $5000$, we set $q=2^{4}=16$, which is close to $\sqrt{1000/4}=15.8$.

To visualize the resulting subsample of Algorithm 1, we generate full data of $2000$ i.i.d. bivariate normal points, truncated in absolute value by 2. The generating distribution has zero mean, unit variance and a correlation of $0.3$ between any two predictors. Figure 2 plots the full data (left), a random subsample (middle), and an IES subsample (right), both subsamples of size $250$. The hyperparameter $q=16$ is used for the IES. Figure 2 clearly shows that predictors in the IES subsample are more uniformly distributed and less correlated than predictors in the random subsample.

6 Simulation Studies

In this section, we evaluate the performance of the IES method through simulation studies. We compare the IES subsample with the random subsample (Rand) and the LowCon method. LowCon is a subsampling method developed in Meng et al. (2020) for smoothing splines and in Meng et al. (2021) for misspecified linear models. It selects a subsample that approximates a prefix space-filling design (Joseph et al., 2015; Lin and Tang, 2015) via nearest neighbor search.

We set the full sample size $N=10000$ and generate values of $p=3$ predictors from two distributional settings:

• Case 1. The predictors follow a truncated multivariate normal $\mathcal{TN}(0,\Sigma,-2,2)$ with mean zero and covariance matrix $\Sigma=(0.3^{\mathbbm{1}{(i\neq j)}}).$ Each predictor lies in $[-2,2]$.
• Case 2. The predictors are generated via a truncated multivariate exponential distribution using the elliptical copula in the R package copula. The covariance matrix $\Sigma$ is the same as in Case 1. The marginal distribution is specified as an exponential with rate one, and is truncated above by $4$ and translated to $[-2,2]$.

The responses are generated by $Y=m(\mathbf{X})+\epsilon$, where

[TABLE]

and $\epsilon$ follows $\mathcal{N}(0,0.25)$.

The effect of model misspecification on IES is also studied, where an additional interaction term $2\ln(4.5+X_{1}X_{2})$ is added to the true regression function in (12) but is not used when training an additive model.

Each setting of predictors is replicated $200$ times, and the three subsampling methods, Rand, LowCon and IES, are performed for each replication with the subsample size $n=1000$. The hyperparameter $q=16$ is used for the IES method. Backfitting algorithm with local linear smoothers is then applied to train an additive model over each subsample. The bandwidth, searched in $\{0.05,0.1,0.15,\dots,0.95\}^{3}$, is chosen via a five-fold cross validation (CV). For the $\hat{m}$ trained over each subsample, we consider two performance measures, namely, the maximum estimation error ${\rm MEE}=\max_{\mathbf{x}\in\mathcal{X}_{test}}|\hat{m}(\mathbf{x})-m(\mathbf{x})|,$ and the average squared error ${\rm ASE}=\sum_{\mathbf{x}\in\mathcal{X}_{test}}\left(\hat{m}(\mathbf{x})-m(\mathbf{x})\right)^{2}/10^{6}.$ The ${\rm MEE}$ is a realization of the maximum risk used in (4) and quantifies the worst performance of $\hat{m}$, and the ${\rm ASE}$ measures the overall performance of $\hat{m}$ over the test domain. The test data $\mathcal{X}_{test}$ are $10^{6}$ grid points with each predictor spanning at 100 evenly spaced points from $-1.8$ to $1.8$.

Figure 3 plots the ${\rm MEE}$ and ${\rm ASE}$ of $\hat{m}$ trained on different subsamples across the 200 replications. The IES consistently allows better estimation of $m$ than the subsamples selected from other methods. Specifically, the MEE plots demonstrate the advantage of the IES in controlling the worst error across the entire domain, and the ASE plots suggest a better overall estimation performance of IES.

Figure 4 depicts the fitted curves of each component function in (12) for each subsample of the full data generated in Case 2. The red curve represents the target centered component function, and the black curve indicates the average fit over the $200$ replications. The grey shaded area is the empirical $95\%$ confidence band. It is clear that the IES method always outperforms random subsampling in allowing a better fit of each component function. When compared with LowCon, the IES performs similarly in terms of average fit, but it performs better in terms of stability (width of the shaded band), especially in the area with low density, i.e. the right tails of all component functions, and when the target function assumes a nonlinear shape, e.g. the turnings areas in the second and third component functions. Figure LABEL:fig:normcompare in Supplementary Materials reveals similar comparison results for the subsamples of the full data generated in Case 1.

The out-performance of IES over LowCon comes from two aspects. Firstly, the IES samples diverse points sequentially and avoids duplicates, while LowCon applies nearest neighborhood search to approximate a prefix space-filling design, which often samples repeatedly on the same observation in the region with scarce data. Duplicated points have bigger weights and increase the modeling instability. Secondly, most space-filling designs target at full dimensional uniformity but may not be uniform when projected to low dimensions. IES targets at one- and two-dimensional uniformity and thus is more suitable for establishing additive models.

Figure 5 shows the boxplots of MEEs and ASEs for the regression function with the misspecified interaction term $2\ln(4.5+X_{1}X_{2})$. The predictors are generated the same as in Case 2. The lower estimation error for IES suggests that its subsamples are less susceptible to model misspecification because of the fact that the predictors in an IES subsample are less dependent. In our particular setting, $X_{3}$ is nearly independent of $X_{1}$ and $X_{2}$ in the IES subsample. Hence, the component function of $X_{3}$ is not affected by the misspecified interaction term of $X_{1}$ and $X_{2}$ and be accurately estimated. The plots of estimated component functions are relegated to Figure LABEL:fig:miscompare in Supplementary Materials to save space.

7 Real Data

We now evaluate the performance of the IES method on the Diamond Price Prediction dataset. The dataset is available from both the R package ggplot2 and https://www.kaggle.com/shivam2503/diamonds. Price along with $9$ predictors of 53,940 diamonds are collected in the data with the goal of building a predictive model for the diamond price. Three discrete quality measures, namely cut, color, and clarity, are dropped, as we focus on continuous predictors. Among continuous predictors, carat, depth (which summarizes information in other left-out predictors) and table, are picked for modeling. The first predictor measures the weight of each diamond and the latter two are specialized shape metrics. Since both carat and price are highly skewed, a log transformation is applied. We train the model

[TABLE]

over selected subsamples via the same backfitting procedure as in Section 6.

7.1 Estimation Performance

Backfitting on the LowCon subsample of this dataset does not converge. Therefore, we only compare the IES with random subsamples. We use the model trained on the full data as a benchmark because the true model is unknown to us. The subsample size is fixed at $n=5000$.

Figure 6 depicts estimated component functions trained on the full sample and subsamples selected by different methods. The span of $x$-axis of each component function reflects its range in the full data. Since a subsample often results in a reduced range of predictors, extrapolation is needed. In this case, we use term-wise nearest neighbor estimation. In Figure 6, the component function of carat has a dominant effect in magnitude with mostly a linear shape. The estimations over an IES subsample and a random subsample are both close to the benchmark, with the IES showing its advantage in the right tail. This confirms that the IES subsample provides better worst-case control in accuracy. The estimation of the other two component functions clearly demonstrates the superiority of IES. The IES effectively captures the information of each component function, even if the function has a complex shape and a relatively weak signal.

Table 1 further compares the performance of IES and random subsamples using measures for estimation and prediction errors. First, same as in Section 6, we calculate MEE and ASE for the regression function over the test data $\mathcal{X}_{test}$, the grid points of size $10^{6}$ that span the range of the full data. The response for $\mathcal{X}_{test}$ is generated using the model trained on the full data. In addition, we calculate the average (AvePredError) and maximum prediction error (MaxPredError) for the observed price in the full data. From Table 1, an IES subsample outperforms a random subsample in minimizing both estimation and prediction errors. An additive model trained on an IES subsample provides more accurate component function estimation and response prediction than the model trained over a random subsample.

7.2 Computation Time

We now report the computational time of IES on the Diamond data. Table 2 lists the computation time of subsampling, CV, and model fitting procedures as well as the total spent time, with their respective standard deviations shown in parenthesis. As shown in Table 2, CV dominates the time consumption for training an additive model, making the modeling on the full data dramatically slow. Training the model on a subsample significantly accelerates the CV and reduces the time to around 8-fold. The IES sampling procedure does take a few more seconds, but this is unimportant compared to the big saving on the time for CV. The total time of IES and Rand are comparable, and it makes sense for IES to be a little slower than Rand to achieve its superior estimation performance.

8 Discussion

We have developed a new subsampling method, called ISE, to accelerate the computation of training an additive model from large data. The ISE selects the subsample that approximates an OA and optimizes the minimax risk of training an additive model by enabling asymptotically independent and uniformly distributed predictors in the selected subsample. Theoretical results have been derived to guarantee the convergence of the backfitting procedure over an ISE subsample for two-dimensional problems. Extensive simulation studies and a real data application demonstrate that ISE outperforms existing subsampling methods in providing accurate estimations of the regression function and each component.

Future works can look into subsampling via OAs with higher strength. The asymptotic property in Theorem 4 can be easily extended to a general number of predictors if the training subsample has a higher strength. In addition, such a subsample achieves higher-order independence among multiple predictors and will allow better estimation of an additive model with interaction terms. Another direction is to consider the performance of IES for a more general family of models, for example, the generalized additive model. We expect that such a subsample will perform well for estimating $g(E[Y])$ for a general link function $g$ because of its independence between predictors and uniform coverage of the data region.

Supplementary Materials

The supplementary materials include the proofs of Theorems 1–4 and additional simulation results.

Appendix A Existence of OA

The existence and construction of OAs have been widely studied in the literature, see, for example, Hedayat et al. (1999) and Dey and Mukerjee (2009) for a comprehensive introduction. Below is a well-known result.

Lemma A1**.**

If $q$ is a prime power and $\lambda$ is a positive integer, then an OA$(\lambda q^{2},p,q,2)$ exists for any $p\leq q+1$.

A construction of OA$(q^{2},q+1,q,2)$ with $q$ being a prime power can be found in Hedayat et al. (1999) (Theorem $3.1$). Stacking $\lambda$ copies of an OA$(q^{2},q+1,q,2)$ provides an OA$(\lambda q^{2},q+1,q,2)$, any $p$ columns of which is an OA$(\lambda q^{2},p,q,2)$.

When $q$ is not a prime power, one may construct OAs from pairwise orthogonal Latin squares. The lemma below comes from this approach.

Lemma A2**.**

Let $q_{1}^{v_{1}}q_{2}^{v_{2}}\cdots q_{u}^{v_{u}}$ be a prime factorization of $q$ and $q_{0}=\min\{q_{i}^{v_{i}}\mid i=1,\dots,u\}$, then an OA$(\lambda q^{2},p,q,2)$ exists for any $p\leq q_{0}+1$.

The result is an immediate consequence of Theorems $8.4$ and $8.28$ in Hedayat et al. (1999). It extends $q$ from prime power to an arbitrary positive integer.

Many other OAs with flexible $p$ and $q$ exist, see http://neilsloane.com/oadir/ for a collection of examples.

Appendix B The unique solution in Theorem 4

Denote the observations in the IES subsample $\{\mathcal{X}^{*}_{opt},\mathbf{y}^{*}_{opt}\}$ by $(x^{*}_{i1},x^{*}_{i2},y^{*}_{i})$ where $i\in\{1,2,\dots,n\}$, $x^{*}_{i1}$ and $x^{*}_{i2}$ are the two predictors, and $y^{*}_{i}$ is the response. Define, for $t=0,1,2$,

[TABLE]

Then define $n\times n$ matrices $\mathcal{S}_{1}=\{[\mathcal{S}_{1}]_{ij}\}_{1\leq i,j\leq n}$ and $\mathcal{S}_{2}=\{[\mathcal{S}_{2}]_{ij}\}_{1\leq i,j\leq n}$ where

[TABLE]

Following Buja et al. (1989) and Opsomer and Ruppert (1997), the bivariate additive model, fitted by local linear smoothers via backfitting algorithm, aims to solve the following estimation equation:

[TABLE]

where $\mathbf{\hat{m}}_{1}=\left(\hat{m}_{1}(x^{*}_{11}),\dots,\hat{m}_{1}(x^{*}_{n1})\right)^{\top}$, $\mathbf{\hat{m}}_{2}=\left(\hat{m}_{2}(x^{*}_{12}),\dots,\hat{m}_{2}(x^{*}_{n2})\right)^{\top}$, $\mathbf{Y}=\left(y^{*}_{1},\dots,y^{*}_{n}\right)^{\top}$, $\mathcal{S}_{1}^{*}=(\mathcal{I}-\mathbf{1}\mathbf{1}^{\top}/n)\mathcal{S}_{1}$, and $\mathcal{S}_{2}^{*}=(\mathcal{I}-\mathbf{1}\mathbf{1}^{\top}/n)\mathcal{S}_{2}$ with $\mathcal{I}$ being the $n\times n$ identity matrix and $\mathbf{1}$ being a $n\times 1$ vector of all ones. The centering constant $\mu$ is estimated separately by $\hat{\mu}=\bar{y}$. The backfitting algorithm on the IES subsample converges to the unique solution

[TABLE]