More efficient approximation of smoothing splines via space-filling basis selection

Consider the nonparametric regression model $y_{i}=\eta(x_{i})+\epsilon_{i}(i=1,\ldots,n)$, where $y_{i}\in\mathbb{R}$ denotes the $i$th observation of the response, $\eta$ represents an unknown function to be estimated, $x_{i}\in\mathbb{R}^{d}$ is the $i$th observation of the predictor variable, and $\{\epsilon_{i}\}_{i=1}^{n}$ are independent and identically distributed random errors with zero mean and unknown variance $\sigma^{2}$. The function $\eta$ can be estimated by minimizing the penalized least squares criterion,

where $J(\eta)$ denotes a quadratic roughness penalty (Wahba, 1990; Wang et al., 2011; Gu, 2013). The smoothing parameter $\lambda$ here administrates the trade-off between the goodness-of-fit of the model and the roughness of the function $\eta$. In this paper, expression (1) is minimized in a reproducing kernel Hilbert space ${\mathcal{H}}$, which leads to a smoothing spline estimate for $\eta$.

Although univariate smoothing splines can be computed in $O(n)$ time (Reinsch, 1967), it takes $O(n^{3})$ time to find the minimizer of (1) when $d\geq 2$. Such a computational cost hinders the use of smoothing splines for large samples. To reduce the computational cost for smoothing splines, extensive efforts have been made to approximate the minimizer of (1) by restricting the estimator $\hat{\eta}$ to a subspace ${\mathcal{H}}_{E}\subset{\mathcal{H}}$. Let the dimension of the space ${\mathcal{H}}_{E}$ be $q$ and the restricted estimator be $\hat{\eta}_{E}$. Compared with the $O(n^{3})$ computational cost of calculating $\hat{\eta}$, the computational cost of $\hat{\eta}_{E}$ can be reduced to $O(nq^{2})$. Along this line of thinking, numerous studies have been developed in recent decades. Luo and Wahba (1997) and Zhang et al. (2004) approximated the minimizer of (1) using variable selection techniques. Pseudosplines (Hastie, 1996) and penalized splines (Ruppert et al., 2009) were also developed to approximate smoothing splines.

Although these methods have already yielded impressive algorithmic benefits, they are usually ad hoc in choosing the value of $q$. The value of $q$ regulates the trade-off between the computational time and the prediction accuracy. One fundamental question is how small $q$ can be in order to ensure the restricted estimator $\hat{\eta}_{E}$ converge to the true function $\eta$ at the identical rate as the full sample estimator $\hat{\eta}$. To answer this question, Gu and Kim (2002); Ma et al. (2015) developed random sampling methods for selecting the basis functions and established the coherent theory for the convergence of the restricted estimator $\hat{\eta}_{E}$. To ensure that $\hat{\eta}_{E}$ has the same convergence rate as $\hat{\eta}$, both methods in Gu and Kim (2002) and Ma et al. (2015) require $q$ be of the order $O\{n^{2/(pr+1)+\delta}\}$, where $\delta$ is an arbitrary small positive number, $p\in[1,2]$ depends on the true function $\eta$, and $r$ depends on the fitted spline. In Gao (1999), it is shown that fewer basis functions are needed to warrant the aforementioned convergence rate if we select the basis functions $\{R(z_{j},\cdot)\}_{j=1}^{q}$, where $\{z_{j}\}_{j=1}^{q}$ are roughly equal-spaced. However, they only provide the theory in the univariate predictor case, and their method cannot be directly applied to multivariate cases.

In this paper, we develop a more efficient computational method to approximate smoothing splines. The distinguishing feature of the proposed method is that it considers the notion of diversity of the selected basis functions. We propose the space-filling basis selection method, which chooses the basis functions with a large diversity, by selecting the ones that correspond to roughly uniformly-distributed observations. When $d\leq pr+1$, we show that the smoothing spline estimator proposed here has the same convergence rate as the full sample estimator, and the order of the essential number of basis function $q$ is reduced to $O\{n^{(1+\delta)/(pr+1)}\}$.

Let ${\mathcal{H}}=\{\eta:J(\eta)<\infty\}$ be a reproducing kernel Hilbert space, where $J(\cdot)$ is a squared semi-norm. Let ${\mathcal{N}}_{J}=\{\eta:J(\eta)=0\}$ be the null space of $J(\eta)$ and assume that ${\mathcal{N}}_{J}$ is a finite-dimensional linear subspace of ${\mathcal{H}}$ with basis $\{\xi_{i}\}_{i=1}^{m}$ in which $m$ is the dimension of ${\mathcal{N}}_{J}$. Let ${\mathcal{H}}_{J}$ be the orthogonal complement of ${\mathcal{N}}_{J}$ in ${\mathcal{H}}$ such that ${\mathcal{H}}={\mathcal{N}}_{J}\oplus{\mathcal{H}}_{J}$. The space ${\mathcal{H}}_{J}$ is a reproducing kernel Hilbert space with $J(\cdot)$ as the squared norm. The reproducing kernel of ${\mathcal{H}}_{J}$ is denoted by $R_{J}(\cdot,\cdot)$. The well-known representer theorem (Wahba, 1990) states that there exist vectors $d=(d_{1},\ldots,d_{m})^{T}\in\mathbb{R}^{m}$ and $c=(c_{1},\ldots,c_{n})^{T}\in\mathbb{R}^{n}$, such that the minimizer of (1) in ${\mathcal{H}}$ is given by $\eta(x)=\sum_{k=1}^{m}d_{k}\xi_{k}(x)+\sum_{i=1}^{n}c_{i}R_{J}(x_{i},x).$ Let $Y=(y_{1},\ldots,y_{n})^{T}$ be the vector of response observations, $S$ be the $n\times m$ matrix with the $(i,j)$th entry $\xi_{j}(x_{i})$, and $R$ be the $n\times n$ matrix with the $(i,j)$th entry $R_{J}(x_{i},x_{j})$. Solving the minimization problem in (1) thus is equivalent to solving

where the smoothing parameter $\lambda$ can be selected based on the generalized cross-validation criterion (Wahba and Craven, 1978). In a general case where $n\gg m$ and $d\geq 2$, the computation cost for calculating $(\widehat{d},\widehat{c})$ in equation (2) is of the order $O(n^{3})$, which is prohibitive when the sample size $n$ is considerable. To reduce this computational burden, one can restrict the full sample estimator $\hat{\eta}$ to a subspace ${\mathcal{H}}_{E}\subset{\mathcal{H}}$, where ${\mathcal{H}}_{E}={\mathcal{N}}_{J}\oplus\mbox{span}\{R_{J}(x_{i}^{*},\cdot),i=1,\ldots,q\}$. Here, ${\mathcal{H}}_{E}$, termed as the effective model space, can be constructed by selecting a subsample $\{x_{i}^{*}\}_{i=1}^{q}$ from $\{x_{i}\}_{i=1}^{n}$. Such an approach is thus called the basis selection method.

Denote a matrix $R_{*}\in\mathbb{R}^{n\times q}$ with the $(i,j)$th entry $R_{J}(x_{i},x_{j}^{*})$ and $R_{**}\in\mathbb{R}^{q\times q}$ with the $(i,j)$th entry $R_{J}(x_{i}^{*},x_{j}^{*})$. The minimizer of (1) in the effective model space ${\mathcal{H}}_{E}$ thus can be written as $\eta_{E}(x)=\sum_{k=1}^{m}d_{k}\xi_{k}(x)+\sum_{i=1}^{q}c_{i}R(x^{*}_{i},x)$ and the coefficients, $d_{E}=(d_{1},\ldots,d_{m})^{T}$ and $c_{E}=(c_{1},\ldots,c_{q})^{T}$ can be obtained through solving

The evaluation of the restricted estimator $\hat{\eta}_{E}$ at sample points thus satisfies $\hat{\eta}_{E}=S\hat{d}_{E}+R_{*}\hat{c}_{E}$, where $\hat{\eta}_{E}=\{\hat{\eta}_{E}(x_{1}),\ldots,\hat{\eta}_{E}(x_{n})\}^{T}$. In a general case where $m\ll q\ll n$, the overall computational cost for the basis selection method is $O(nq^{2})$, which is a significant reduction compared with $O(n^{3})$. Recall that the value of $q$ controls the trade-off between the computational time and the prediction accuracy. To ensure that $\hat{\eta}_{E}$ converges to the true function $\eta$ at the same rate as $\hat{\eta}$, researchers showed that the essential number of basis functions $q$ needs to be of the order $O\{n^{2/(pr+1)+\delta}\}$, where $\delta$ is an arbitrary small positive number (Kim and Gu, 2004; Ma et al., 2015). In the next section, we present the proposed space-filling basis selection method, which reduces such an order to $O\{n^{(1+\delta)/(pr+1)}\}$.

Recall that the data points are assumed to be generated from the uniform distribution on $[0,1]^{d}$. Thus, for each coordinate $x$, the corresponding marginal density $f_{X}(\cdot)=1$. We define that $V(g)=\int_{[0,1]^{d}}g^{2}\mbox{d}x$. The following four regularity conditions are required for the asymptotic analysis, and the first three are the identical conditions employed by Ma et al. (2015), in which one can find more technical discussions.

Condition 1. The function $V$ is completely continuous with respect to $J$;

Condition 2. for some $\beta>0$ and $r>1$, $\rho_{\nu}>\beta\nu^{r}$ for sufficiently large $\nu$;

Condition 3. for all $\mu$ and $\nu$, $\mbox{var}\{\phi_{\nu}(x)\phi_{\mu}(x)\}\leq C_{1}$, where $\phi_{\nu}$, $\phi_{\mu}$ are the eigenfunctions associated with $V$ and $J$ in ${\mathcal{H}}$, $C_{1}$ denotes a positive constant;

Condition 4. for all $\mu$ and $\nu$, $\mathcal{V}(g_{\nu,\mu})\leq C_{2}$, where $\mathcal{V}(\cdot)$ denotes the total variation, $g_{\nu,\mu}(x)=\phi_{\nu}(x)\phi_{\mu}(x)$, and $C_{2}$ represents a positive constant. The total variation here is defined in the sense of Hardy and Krause (Owen, 2003). As a specific case when $d=1$, the total variation $\mathcal{V}(g)=\int|g^{\prime}(x)|\mbox{d}x$, revealing that a smooth function displays a small total variation. Intuitively, the total variation measures how wiggly the function $g$ is.

Condition 1 indicates that there exist a sequence of eigenfunctions $\phi_{\nu}\in{\mathcal{H}}$ and the associated sequence of eigenvalues $\rho_{\nu}\uparrow\infty$ satisfying $V(\phi_{\nu},\phi_{\mu})=\delta_{\nu\mu}$ and $J(\phi_{\nu},\phi_{\mu})=\rho_{\nu}\delta_{\nu\mu}$, where $\delta_{\nu\mu}$ is the Kronecker delta. The growth rate of the eigenvalues $\rho_{\nu}$ dictates how fast $\lambda$ should approach to 0, and further what the convergence rate of smoothing spline estimates is (Gu, 2013). Notice that the eigenfunctions $\phi_{\nu}$s have a close relationship with the Demmler-Reinsch basis, which are orthogonal vectors representing $l_{2}$ norm  (Ruppert, 2002). The eigenfunctions $\phi_{\nu}$s can be calculated explicitly in several specific scenarios. For instance, $\phi_{\nu}$s are the sine and cosine functions when $J(\eta)=\int_{0}^{1}(\eta^{{}^{\prime\prime}})^{2}\mbox{d}x$, where $\eta$ denotes a periodic function on $[0,1]$. For more details on the construction of $\phi_{\nu}$s can be found in Section 9.1 of Gu (2013).

Combining Lemma 1 and Theorem 1 and set $h=g_{\nu,\mu}$, $\mathcal{T}_{q}={\mathcal{X}}_{q}^{*}$ yields the following lemma.

Lemma 2 shows the advantage of $\{x_{i}^{*}\}_{i=1}^{q}$, the subsample selected by the proposed method, over a randomly selected subsample $\{x_{i}^{+}\}_{i=1}^{q}$. To be specific, as a direct consequence of Condition 3, we have $E[\int_{[0,1]^{d}}\phi_{\nu}\phi_{\mu}\mbox{d}x-\sum_{j=1}^{q}\phi_{\nu}(x_{j}^{+})\phi_{\mu}(x_{j}^{+})/q]^{2}=O(q^{-1}),$ for all $\mu$ and $\nu$. Lemma 2 therefore implies the subsample ${\mathcal{X}}_{q}^{*}$ can more efficiently approximate the integration $\int_{[0,1]^{d}}\phi_{\nu}\phi_{\mu}\mbox{d}x$, for all $\mu$ and $\nu$. We now present our main theoretical result.

It is shown in Theorem 9.17 of Gu (2013) that the full sample smoothing spline estimator $\hat{\eta}$ has the convergence rate, $(V+\lambda J)(\hat{\eta}-\eta_{0})=O_{p}(n^{-pr/(pr+1)})$ under some regularity conditions. Theorem 2 thus states that proposed estimator $\hat{\eta}_{E}$ achieves the same convergence rate as the full sample estimator $\hat{\eta}$, under two extra conditions imposed on $q$ (a) $q=O(n^{1/d})$, and (b) $q^{1-\delta}\lambda^{1/r}\rightarrow\infty$ as $\lambda\rightarrow 0$. Moreover, Theorem 2 indicates that to achieve the identical convergence rate as the full sample estimator $\hat{\eta}$, the proposed approach requires a much smaller number of basis functions, in the case when $\lambda\asymp n^{-r/(pr+1)}$. $q^{1-\delta}\lambda^{1/r}\rightarrow\infty$ indicates an essential choice of $q$ for the proposed estimator should satisfy $q=O\{n^{(1+\delta)/(pr+1)}\}$, when $\lambda\asymp n^{-r/(pr+1)}$. For comparison, both the random basis selection (Gu and Kim, 2002) and the adaptive basis selection method (Ma et al., 2015) require the essential number of basis functions to be $q=O\{n^{2/(pr+1)+\delta}\}$. As a result, the proposed estimator is more efficient since it reduces the order of the essential number of basis functions.

Given $q=O(n^{1/d})$, when $d\leq pr+1$, it follows naturally when $q^{1-\delta}\lambda^{1/r}\rightarrow\infty$ is satisfied. Otherwise, when $d>pr+1$, $q=O(n^{1/d})$ becomes sufficient but not necessary for $q^{1-\delta}\lambda^{1/r}\rightarrow\infty$. We thus stress that the essential number of basis functions for the proposed method, $q=O\{n^{(1+\delta)/(pr+1)}\}$, can only be achieved when $d\leq pr+1$. The parameter $p$ in Theorem 2 is closely associated with the true function $\eta_{0}$ and will affect the convergence rate of the proposed estimator. Intuitively, the larger the $p$ is, the smoother the function $\eta_{0}$ will be. For $p\in[1,2]$, the optimal convergence rate of $(V+\lambda J)(\hat{\eta}_{E}-\eta_{0})$ falls in the interval $[O_{p}(n^{-r/(r+1)}),O_{p}(n^{-2r/(2r+1)})]$. To the best of our knowledge, the problem of selecting the optimal $p$ has rarely been studied, and one exception is Serra and Krivobokova (2017), where the author studied such a problem in one-dimensional cases. Serra and Krivobokova (2017) provided a Bayesian approach for selecting the optimal parameter, named $\beta$, which is known to be proportional to $p$. Nevertheless, since the constant $\beta/p$ is usually unknown, such an approach still cannot be used to select the optimal $p$ in practice. Furthermore, whether such an approach can be extended to the high-dimensional cases remains unclear.

For the dimension of the effective model space $q$, a suitable choice is $q=n^{(1+\delta)/(4p+1)+\delta}$ in the following two cases. Case 1. Univariate cubic smoothing splines with the penalty $J(\eta)=\int_{0}^{1}(\eta^{\prime\prime})^{2}$, $r=4$ and $\lambda\asymp n^{-4/(4p+1)}$; Case 2. Tensor-product splines with $r=4-\delta^{*}$, where $\delta^{*}>0$. For $p\in[1,2]$, the dimension roughly lies in the interval $(O(n^{1/9}),O(n^{1/5}))$.

To assess the performance of the proposed space-filling basis selection method, we carry out extensive analysis on simulated datasets. We report some of them in this section. We compare the proposed method with uniform basis selection and adaptive basis selection, and report both prediction error and running time.

The following four functions on $[0,1]$ (Lin and Zhang, 2006) are used as the building blocks in our simulation study, $g_{1}(t)=t$, $g_{2}(t)=(2t-1)^{2}$, $g_{3}(t)=\sin(2\pi t)/\{2-\sin(2\pi t)\}$, and $g_{4}(t)=0.1\sin(2\pi t)+0.2\cos(2\pi t)+0.3\sin(2\pi t)^{2}+0.4\cos(2\pi t)^{3}+0.5\sin(2\pi t)^{3}$. In addition, we also use the following functions on $[0,1]^{2}$ (Wood, 2003) as the building blocks,

where $\sigma_{1}=0.3$ and $\sigma_{2}=0.4$. The signal-to-noise ratio (SNR), defined as $\mbox{var}\{\eta(X)\}/\sigma^{2}$, is set to be at two levels: 5 and 2. We generate replicated samples with sample sizes $n=\{2^{10},2^{11},\ldots,2^{14}\}$ and dimensions $d=\{2,4,6\}$ uniformly on $[0,1]^{p}$ from the following four regression settings,

(1) A 2-d function $g_{1}(x_{1}x_{2})+g_{2}(x_{2})+g_{3}(x_{1})+g_{4}(x_{2})+g_{3}\{(x_{1}+x_{2})/2\}$;

(2) A 2-d function $h_{1}(x_{1},x_{2})+h_{2}(x_{1},x_{2})$;

(3) A 4-d function $g_{1}(x_{1})+g_{2}(x_{2})+g_{3}(x_{3})+2g_{1}\{(x_{1}+x_{4})/2\}+2g_{2}\{(x_{2}+x_{3})/2\}+2g_{3}\{(x_{1}+x_{3})/2\}$;

(4) A 6-d function $h(x_{1},x_{2})+h(x_{1},x_{5})$.

In the simulation, $q$ is set to be $5n^{2/9}$ and $10n^{1/9}$, based on the asymptotic results. To combat the curse of dimensionality, we fit smoothing spline analysis of variance models with all main effects and two-way interactions. The prediction error is measured by the mean squared error (MSE), defined as $[\sum_{i=1}^{n_{0}}\{\hat{\eta}_{E}(t_{i})-\eta_{0}(t_{i})\}^{2}]/n_{0}$, where $\{t_{i}\}_{i=1}^{n_{0}}$ denotes an independent testing dataset uniformly generated on $[0,1]^{p}$ with $n_{0}=5000$. The max projection design (Joseph et al., 2015) is used to generate design points in Step 1 of the proposed method. Our empirical studies suggest that the Sobol sequence and other space-filling techniques, e.g., the Latin hypercube design (Pukelsheim, 2006) and the uniform design (Fang et al., 2000), also yield similar performance.

Figure 3 shows the MSE against the sample size on the log-log scale. Each column presents the results of a function setting as described above. We set the signal-to-noise ratio to be five and two in the upper row and the lower row, respectively. We use solid lines for the proposed method, dotted lines for adaptive basis selection method, and dashed lines for uniform basis selection method. The number of the basis functions $q$ is $5n^{2/9}$ and $10n^{1/9}$ for the lines with triangles and the lines with circles, respectively. The vertical bars represent standard error bars obtained from 20 replicates. The full sample estimator is omitted here due to the high computation cost. It is observed that the space-filling basis selection method provides more accurate smoothing spline predictions than the other two methods in almost all settings. Notably, the lines with circles for the space-filling basis selection method displays a linear trend similar to the lines with triangles for the other two methods. This observation reveals the proposed estimator yields a faster convergence rate than the other two methods.

More simulation results can be found in the Supplementary Material, in which we consider the regression functions that exhibit several sharp peaks. In those cases, the results suggest that both the space-filling basis selection method and the adaptive basis selection method outperform the uniform basis selection, whereas neither the space-filling basis selection method nor the adaptive basis selection method dominates each other. Moreover, the proposed space-filling basis selection method outperforms the adaptive basis selection method as the sample size $n$ gets larger.

Table 1 summarizes the computing times of model-fitting using all methods on a synthetic dataset with $n=2^{14}$ and $q=5n^{2/9}$. The simulation is replicated for 20 runs using a computer with an Intel 2.6 GHz processor. In Table 1, UNIF, ABS, and SBS represent the uniform basis selection method, the adaptive basis selection method, and the proposed space-filling basis selection method, respectively. The time for calculating the smoothing parameter is not included. The result for the full basis smoothing spline estimator is omitted here due to the huge computational cost. The computational time for generating a set of design points, i.e., Step 1 in the proposed algorithm, is not included since the design points can be generated beforehand. It is observed that the computing time of the proposed method is comparable with that of the other two basis selection methods under all settings. Combining such an observation with the result in Fig. 3, it is concluded that the proposed method can achieve a more accurate prediction, without requiring much more computational time.

The problem of measuring total column ozone has attracted significant attention for decades. Ozone depletion facilitates the transmission of ultraviolet radiation (290–400 nm wavelength) through the atmosphere and causes severe damage to DNA and cellular proteins that are involved in biochemical processes, affecting growth and reproduction. Statistical analysis of total column ozone data has three steps. In the first step, the raw satellite data (level 1) are retrieved by NASA. Subsequently, NASA calibrates and preprocesses the data to generate spatially and temporally irregular total column ozone measurements (level 2). Finally, the level 2 data are processed to yield the level 3 data, which are the daily and spatially regular data product released extensively to the public.

We fit the nonparametric model $y_{ij}=\eta(x_{\langle 1\rangle i},x_{\langle 2\rangle j})+\epsilon_{ij}$ to a level 2 total column ozone dataset ($n$=173,405) compiled by Cressie and Johannesson (2008). Here, $y_{ij}$ is the level 2 total column ozone measurement at the $i$th longitude, i.e., $x_{\langle 1\rangle i}$, and the $j$th latitude, i.e., $x_{\langle 2\rangle i}$, and $\epsilon_{ij}$ represent the independent and identically distributed random errors. The heat map of the raw data is presented in the Supplementary Material. The thin-plate smoothing spline is used for the model-fitting, and the proposed method is employed to facilitate the estimation. The number of basis functions is set to $q=20n^{2/9}\approx 292$. The design points employed in the proposed basis selection method are yielded from a Sobol sequence (Dutang and Savicky, 2013). The heat map of the predicted image on a $1^{\circ}\times 1^{\circ}$ regular grid is presented in Fig. 4. It is seen that the total column ozone value decreases dramatically to form the ozone hole over the South Pole, around –$55^{\circ}$ latitudinal zone.

We now report computing times of the model-fitting that are performed on the identical laptop computer for the simulation studies. The computational times, in CPU seconds, are presented in parentheses, including basis selection (0.1s), model fitting (129s), and prediction (21s). Further comparison between the proposed method and other basis selection methods on this dataset can be found in the Supplementary Material.

This study was partially supported by the National Science Foundation and the National Institutes of Health.