Consider the nonparametric regression model y i = η(x i ) + ǫ i (i = 1, . . . , n), where y i ∈ R denotes the ith observation of the response, η represents an unknown function to be estimated, x i ∈ R d is the ith observation of the predictor variable, and {ǫ i } n i=1 are independent and identically distributed random errors with zero mean and unknown variance σ 2 . The function η can be estimated by minimizing the penalized least squares criterion,
where J(η) denotes a quadratic roughness penalty [Wahba, 1990 , Wang et al., 2011 , Gu, 2013 . The smoothing parameter λ here administrates the trade-off between the goodness-of-fit of the model and the roughness of the function η. In this paper, expression (1) is minimized in a reproducing kernel Hilbert space H, which leads to a smoothing spline estimate for η.
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 ≥ 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η to a subspace H E ⊂ H. Let the dimension of the space H E be q and the restricted estimator beη E . Compared with the O(n 3 ) computational cost of calculatingη, the computational cost ofη 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η E converge to the true function η at the identical rate as the full sample estimatorη. 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η E . To ensure thatη E has the same convergence rate asη, both methods in Gu and Kim [2002] and Ma et al. [2015] require q be of the order O{n 2/(pr+1)+δ }, where δ is an arbitrary small positive number, p ∈ [1, 2] depends on the true function η, 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 , ·)} q j=1 , where {z j } q j=1 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 spacefilling 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 ≤ 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+δ)/(pr+1) }.
The space H J is a reproducing kernel Hilbert space with J(·) as the squared norm. The reproducing kernel of H J is denoted by R J (·, ·). The well-known representer theorem [Wahba, 1990] states that there exist vectors d = (d 1 , . . . , d m ) T ∈ R m and c = (c 1 , . . . , c n ) T ∈ R n , such that the minimizer of (1) in H is given by
. Let Y = (y 1 , . . . , y n ) T be the vector of response observations, S be the n × m matrix with the (i, j)th entry ξ j (x i ), and R be the n × 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 λ can be selected based on the generalized cross-validation criterion [Wahba and Craven, 1978] . In a general case where n ≫ m and d ≥ 2, the computation cost for calculating ( d, 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η to a subspace
Here, H E , termed as the effective model space, can be constructed by selecting a subsample
Such an approach is thus called the basis selection method.
Denote a matrix R * ∈ R n×q with the (i, j)th entry R J (x i , x * j ) and R * * ∈ R q×q with the (i, j)th entry R J (x * i , x * j ). The minimizer of (1) in the effective model space H E thus can be written as
i , x) and the coefficients, d E = (d 1 , . . . , d m ) T and c E = (c 1 , . . . , c q ) T can be obtained through solving
The evaluation of the restricted estimatorη E at sample points thus satisfiesη E = Sd E +R * ĉE , whereη E = {η E (x 1 ), . . . ,η E (x n )} T . In a general case where m ≪ q ≪ 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η E converges to the true function η at the same rate asη, researchers showed that the essential number of basis functions q needs to be of the order O{n 2/(pr+1)+δ }, where δ is an arbitrary small positive number 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+δ)/(pr+1) }.
3 Space-filling basis selection
To motivate the development of the proposed method, we first re-examine the ensemble learning methods, which are wellknown in statistics and machine learning [Dietterich, 2002 , Rokach, 2010 . To achieve better predictive performance than a single learner 1 , ensemble learning methods first build a committee which consists of a number of different learners, then aggregate the predictions of these learners in the committee. The aggregation is usually achieved by employing the majority vote or by calculating a weighted average. The diversity among the learners in the committee holds the key to the success of the ensemble learning methods. A large diversity in the committee yields a better performance of ensemble learning methods [Kuncheva and Whitaker, 2003 ].
The restricted smoothing spline estimatorη E can be considered as an ensemble learning method. In particular, the prediction ofη E is conducted by taking a weighted average of the predictions of the selected basis functions R J (x * i , ·), i ∈ {1, . . . , q}, in addition to the basis functions in the null space N J . Inspired by Kuncheva and Whitaker [2003] , we propose to select a subsample {x * i } q i=1 , such that the diversity among the basis functions {R J (x * i , ·)} q i=1 is as large as possible. One crucial question is how to measure the diversity among a set of basis functions. Notice that adjacent data points, i.e.,
. These observations inspire us to select a set of basis functions
are as uniformly-distributed as possible. The uniformly-distributed property, usually termed as the space-filling property in the experimental design literature [Pukelsheim, 2006] , can be systematically measured by the star discrepancy.
Since the star discrepancy is defined for data in [0, 1] d , in practice, we need to map the data with arbitrary distribution to this domain. Suppose X n = {x i } n i=1 are independent and identically distributed observations generated from the cumulative distribution function F with bounded support D ⊂ R d . Suppose τ is a transformation, such that {τ (x i )} n i=1 has the uniform distribution on [0, 1] d . In a simple case where d = 1 and F is known, we can find the transformation τ by setting τ = F . In a more general case where d > 1 and F is unknown, the transformation τ can be calculated via the optimal transport theory [Villani, 2008] . However, finding the exact solution via the optimal transport theory can be time-consuming. Instead, one may approximate the transformation τ using the iterative transformation approach [Pukelsheim, 2006] or the sliced optimal transport map approach [Rabin et al., 2011] . The computational cost of these two approaches is of the order O{Kn log(n)}, where K denotes the number of iterations [Kolouri et al., 2018 , Cuturi and Doucet, 2014 , Bonneel et al., 2015 . Such a computational cost is negligible compared with the computational cost of the proposed method. In practice, the data can always be preprocessed using the transformation of τ . Without loss of generality, X n may be assumed to be independent and identically distributed observations generated from the uniform distribution on [0, 1] d .
be a hyper-rectangle and 1{·} be the indicator function. The local discrepancy and the star discrepancy are defined as follows [Fang et al., 2005 , Pukelsheim, 2006 .
The star discrepancy corresponding to X q is defined as D * (X q ) = sup a∈[0,1] d D(X q , a).
The local discrepancy D(X q , a) measures the difference between the volume of the hyper-rectangle [0, a) and the fraction of the data points located in [0, a). The local discrepancy is illustrated in the left panel of Fig. 1 . The star discrepancy D * (X q ) calculates the supreme of all the local discrepancy over a ∈ [0, 1] d . In other words, the smaller the D * (X q ) is, the more space-filling the data points X q are [Fang et al., 2005] . Chung [1949] showed that when X q is generated from the uniform distribution in [0, 1] d , D * (X q ) converges to 0 with the order of convergence O[{log log(q)/q} 1/2 ]. Faster convergence rate can be achieved using space-filling design methods [Pukelsheim, 2006] and the low-discrepancy sequence method [Halton, 1960 , Sobol, 1967 , Owen, 2003 ]. The space-filling design methods, developed in the experimental design literature, aim to generate a set of design points that can cover the space as uniformly as possible. For further details, please refer to Wu and Hamada [2011], Fang et al. [2005] , Pukelsheim [2006] . The lowdiscrepancy sequence method Such a method is frequently applied in the field of quasi-Monte Carlo and is extensively employed for numerical integration. For a Sobol sequence S q , one type of low-discrepancy sequence, it is known that D * (S q ) is of the order O{log(q) d /q}, which is roughly the square order D * (X q ) for fixed d. For more in-depth discussions on the quasi-Monte Carlo methods, see e.g., Lemieux [2009] , Leobacher and Pillichshammer [2014] , Dick et al. [2013] or Chapter 5 in Glasserman [2013] and references therein.
Existing studies suggested that space-filling property can be exploited to achieve a fast convergence rate for numerical integration and response surface estimation [Fang et al., 2005 , Pukelsheim, 2006 . These results inspire us to select the space-filling basis functions in smoothing splines. Unfortunately, the existing techniques of space-filling design cannot be applied to our basis selection problem directly due to the following fact. The design space in space-filling design methods is usually continuous, whereas our sample space {x i } n i=1 is finite and discrete. We propose an algorithm to overcome the barrier.
We shall develop a space-filling basis selection method, in which we select the space-filling data points in a computationallyattractive manner. First, a set of design points S q = {s i } q i=1 ∈ [0, 1] d are generated, either using low-discrepancy sequence or space-filling design methods. Subsequently, the nearest neighbor x * i is selected for each s i , from the sample points
can approximate the design points S q well, if x * i is sufficiently close to s i , for i = 1, . . . q. The proposed method is summarized as follows.
• Step 1. Generate a set of design points
Step 2. Select the nearest neighbor for each design point s i from {x i } n i=1 . Let the selected data points be {x * i } q i=1 ,.
•
Step 3. Minimize the penalized least squares criterion (1) over the following effective model space
The proposed algorithm is illustrated through a toy example in the right panel of Fig. 1 . One hundred data points (gray dots) are generated from the uniform distribution in [0, 1] 2 , and a set of design points (black triangles) are generated through the max projection design [Joseph et al., 2015] , a recently developed space-filling design method. The nearest neighbor to each design point is selected (circle). It is observed that the selected subsample is space-filling since it can effectively approximate the design points. In Fig. 2 , the proposed space-filling basis selection method is compared with the uniform basis selection method [Gu and Kim, 2002] and the adaptive basis selection method [Ma et al., 2015] using a two-dimensional toy example. We generate 5000 data points from the uniform distribution in [0, 1] 2 . The leftmost panel in Fig. 2 presents the heat map for the true response surface y = sin{5(x 1 + x 2 )}. The dimension of the effective model space q is set to be 5 × (5000) 2/9 ≈ 33, for all basis selection methods. The selected basis functions are labeled as solid dots in each panel. The right three panels of Fig. 2 plot the heat maps of the spline estimates of all three basis selection methods. In the uniform basis selection method, the default random number generator in R is employed to select the basis functions. It is observed that the selected points are not uniformly distributed. This is a very common phenomenon for uniform basis selection since the randomly selected points do not necessarily look uniformly-distributed, especially when the number of selected points is small. In contrast, it is observed that the basis functions selected by the proposed method are space-filling. Using the space-filling design techniques, the proposed method overcomes the pitfall of uniform basis selection method and uniformly-distribute the selected points. The true response can be better estimated using the proposed method than using other methods. Now we calculate the computational cost of the proposed method. In Step 1, the design points can be generated beforehand; thus, the computational cost in Step 1 can be ignored. For the nearest neighbor search in Step 2, we employ the k-d tree method, which takes O{n log(n)} flops [Bentley, 1975, Wald and Havran, 2006 ]. The computational cost of this step can be further reduced if we are willing to sacrifice the accuracy of the searching results, e.g., using those approximate nearest neighbor searching algorithms [Altman, 1992 , Arya et al., 1994 . For Step 3, computing the smoothing spline estimates in the restricted space H E is of the order O(nq 2 ), as discussed in Section 2.2. In summary, the overall computational cost for the space-filling basis selection method is of the order O(nq 2 ).
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 (·) = 1. We define that V (g) = [0,1] d g 2 dx. 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 2. for some β > 0 and r > 1, ρ ν > βν r for sufficiently large ν;
Condition 3. for all µ and ν, var{φ ν (x)φ µ (x)} ≤ C 1 , where φ ν , φ µ are the eigenfunctions associated with V and J in H, C 1 denotes a positive constant;
Condition 4. for all µ and ν, V(g ν,µ ) ≤ C 2 , where V(·) denotes the total variation, g ν,µ (x) = φ ν (x)φ µ (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 V(g) = |g ′ (x)|dx, 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 φ ν ∈ H and the associated sequence of eigenvalues ρ ν ↑ ∞ satisfying V (φ ν , φ µ ) = δ νµ and J(φ ν , φ µ ) = ρ ν δ νµ , where δ νµ is the Kronecker delta. The growth rate of the eigenvalues ρ ν dictates how fast λ should approach to 0, and further what the convergence rate of smoothing spline estimates is [Gu, 2013] . Notice that the eigenfunctions φ ν s have a close relationship with the Demmler-Reinsch basis, which are orthogonal vectors representing l 2 norm [Ruppert, 2002] . The eigenfunctions φ ν s can be calculated explicitly in several specific scenarios. For instance, φ ν s are the sine and cosine functions when J(η) = 1 0 (η ′′ ) 2 dx, where η denotes a periodic function on [0, 1]. For more details on the construction of φ ν s can be found in Section 9.1 of Gu [2013] .
We now present our main theoretical results, and all the proofs are relegated to the Supplementary Material. For a set of design points S q of size q, we now assume the star discrepancy D * (S q ) converges to zero at the rate of O{log(q) d /q}, or O{q −(1−δ) } for an arbitrary small positive number δ. Such a convergence rate is warranted if S q is generated from a lowdiscrepancy sequence or space-filling design methods, as discussed in Section 3.1. Recall that the proposed method aims to select a subsample that is space-filling, and the key to success depends on whether the chosen subsample X * q can effectively approximate the design points S q . The following lemma bounds the difference between X * q and S q in terms of star discrepancy.
Lemma 1 states that the selected subsample X * q can effectively approximate the design points S q in the sense that the convergence rate of D * (X * q ) is similar to that of D * (S q ). The following theorem is the Koksma-Hlawka inequality, which will be used in proving our main theorem. See Kuipers and Niederreiter [2012] for a proof.
Theorem 1 (Koksma-Hlawka inequality) Let T q = {t 1 , . . . , t q } be a set of data points in [0, 1] d , and h be a function defined on [0, 1] d with bounded total variation V(h). We have
Combining Lemma 1 and Theorem 1 and set h = g ν,µ , T q = X * q yields the following lemma.
Lemma 2 If q = O(n 1/d ), under Condition 4, for all µ and ν, we have
Lemma 2 shows the advantage of {x * i } q i=1 , the subsample selected by the proposed method, over a randomly selected subsample
To be specific, as a direct consequence of Condition 3, we have E
, for all µ and ν. Lemma 2 therefore implies the subsample X * q can more efficiently approximate the integration [0,1] d φ ν φ µ dx, for all µ and ν. We now present our main theoretical result.
Theorem 2 Suppose i ρ p i V (η 0 , φ i ) 2 < ∞ for some p ∈ [1, 2] and δ is an arbitrary small positive number. Under Conditions 1 -4, q = O(n 1/d ), as λ → 0 and q 1−δ λ 1/r → ∞, we have (V + λJ)(η E − η 0 ) = O p (n −1 λ −1/r + λ p ). In particular, if λ ≍ n −r/(pr+1) , the estimator achieves the optimal convergence rate (V + λJ)(η E − η 0 ) = O p {n −pr/(pr+1) }.
It is shown in Theorem 9.17 of Gu [2013] that the full sample smoothing spline estimatorη has the convergence rate, (V + λJ)(η − η 0 ) = O p (n −pr/(pr+1) ) under some regularity conditions. Theorem 2 thus states that proposed estimatorη E achieves the same convergence rate as the full sample estimatorη, under two extra conditions imposed on q (a) q = O(n 1/d ), and (b) q 1−δ λ 1/r → ∞ as λ → 0. Moreover, Theorem 2 indicates that to achieve the identical convergence rate as the full sample estimatorη, the proposed approach requires a much smaller number of basis functions, in the case when λ ≍ n −r/(pr+1) . q 1−δ λ 1/r → ∞ indicates an essential choice of q for the proposed estimator should satisfy q = O{n (1+δ)/(pr+1) }, when λ ≍ 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)+δ }. 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 ≤ pr + 1, it follows naturally when q 1−δ λ 1/r → ∞ is satisfied. Otherwise, when d > pr + 1, q = O(n 1/d ) becomes sufficient but not necessary for q 1−δ λ 1/r → ∞. We thus stress that the essential number of basis functions for the proposed method, q = O{n (1+δ)/(pr+1) }, can only be achieved when d ≤ pr+1. The parameter p in Theorem 2 is closely associated with the true function η 0 and will affect the convergence rate of the proposed estimator. Intuitively, the larger the p is, the smoother the function η 0 will be. For p ∈ [1, 2], the optimal convergence rate of (V + λJ)(η E − η 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 β, which is known to be proportional to p. Nevertheless, since the constant β/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+δ)/(4p+1)+δ in the following two cases. Case 1. Univariate cubic smoothing splines with the penalty J(η) = 1 0 (η ′′ ) 2 , r = 4 and λ ≍ n −4/(4p+1) ; Case 2. Tensor-product splines with r = 4 − δ * , where δ * > 0. For p ∈ [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πt)/{2−sin(2πt)}, and g 4 (t) = 0.1 sin(2πt)+0.2 cos(2πt)+0.3 sin(2πt) 2 +0.4 cos(2πt) 3 + 0.5 sin(2πt) 3 . In addition, we also use the following functions on [0, 1] 2 [Wood, 2003] as the building blocks,
where σ 1 = 0.3 and σ 2 = 0.4. The signal-to-noise ratio (SNR), defined as var{η(X)}/σ 2 , is set to be at two levels: 5 and 2. We generate replicated samples with sample sizes n = {2 10 , 2 11 , . . . , 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
(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 [
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 : Simulation under different settings (from left to right) with SNR being five (the upper row) and two (the lower row). The prediction errors are plotted versus sample size n The vertical bars are standard error bars obtained from 20 replicates. The solid lines, dotted lines, and dashed lines refer to the space-filling basis selection (SBS), the adaptive basis selection (ABS), and the uniform basis selection (UNIF), respectively. The lines with triangles and those with circles represent q = 5n 2/9 and q = 10n 1/9 , respectively. 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 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 = η(x 1 i , x 2 j ) + ǫ 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 ith longitude, i.e., x 1 i , and the jth latitude, i.e., x 2 i , and ǫ 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 ≈ 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 • × 1 • 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 • 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.