nonparametric functional estimation, RKHS, deep learning, computationally efficient learning, hyperparameter scaling, convergence

1 Introduction↩︎

Machine learning models have recently seen a rapid increase in complexity, particularly with the success of large-scale deep neural networks across language, vision, multimodal learning, and generative modeling [1]–[6]. Despite their effectiveness, developing such models often requires substantial computational resources and human effort, especially in the presence of large data volumes, high dimensionality, and complex data structures [3], [7], [8]. As models and datasets continue to grow, it becomes increasingly important to develop principled strategies that use inexpensive preliminary models to guide the development of more complex models.

Traditionally, knowledge distillation (KD) has served as a widely used framework for mitigating computational costs [9]. In the standard KD paradigm, a large and well-trained teacher model transfers information to a smaller student model, allowing the student to retain much of the teacher’s predictive performance with lower computational cost [10]–[15]. While this teacher-to-student paradigm has been highly effective for model compression and efficient deployment, it leaves a critical bottleneck unresolved: standard KD relies on the prerequisite of a fully trained, highly capable teacher model, offering no relief for the exorbitant initial cost of developing the teacher itself.

To address this fundamental limitation in the training phase, we ask whether the mechanics of knowledge transfer can be inverted. Specifically, can a small and computationally efficient student model provide useful information for training a larger teacher model? We propose Knowledge Cascade (KCas), a reverse knowledge distillation framework in which information learned from a student model is transferred upward to guide the teacher model. In KCas, the student is not intended to replace the teacher; rather, it serves as a low-cost preliminary model that extracts useful information for the teacher at a substantially reduced computational cost. While transferring knowledge upward to a model with strictly greater representational capacity may initially seem counterintuitive, the key observation behind KCas is that the student does not need to learn everything the teacher will eventually learn. It only needs to extract transferable information that remains structurally useful when moving from a simpler setting to a more complex one. This dynamic is analogous to a statistical pilot study, where a small preliminary experiment cannot replace the full study but provides valuable guidance for conducting it more efficiently. KCas exploits this principle by learning from a low-cost student model and transferring the insights to the teacher through statistical theory or empirical scaling laws.

We formalize this principle first within the context of nonparametric multivariate functional estimation in reproducing kernel Hilbert spaces (RKHS). In this setting, a subsample-based estimator naturally serves as the student and the full-sample estimator as the teacher, since the effective complexity of the estimator grows with sample size [16]. The main computational bottleneck is selecting multiple smoothing parameters, whose number and search cost increase rapidly with the number of predictors and interaction terms. KCas addresses this bottleneck by selecting smoothing parameters on a small subsample and transferring them to the full sample through asymptotic sample-size scaling. The teacher model is then fitted on the full data using these transferred parameters, avoiding costly full-sample tuning while retaining the statistical benefit of using all observations. Our theory and experiments show that KCas substantially reduces computation and can sometimes improve estimation accuracy by avoiding unstable or overly adaptive full-sample tuning.

While smoothing spline ANOVA models provide the primary theoretical anchor for this paper, they represent just one concrete instance of a fundamental principle: whenever an appropriate scaling relationship exists, information learned from a low-cost student model can be transferred to guide a more expensive teacher model. To demonstrate this versatility, we extend the KCas to kernel density estimation, where bandwidths follow classical asymptotic scaling laws, and to deep learning, where hyperparameters selected on compact student networks can be used to guide larger teacher networks. Ultimately, these extensions position KCas not merely as a specialized technique for RKHS, but as a general, theoretically grounded paradigm for reverse knowledge transfer across diverse machine learning domains.

Our contributions:

1. We introduce Knowledge Cascade (KCas), a reverse knowledge distillation framework in which small, computationally efficient student models guide larger and more complex teacher models, providing a new perspective beyond conventional teacher-to-student distillation.

2. We show that student-to-teacher knowledge transfer can be principled when supported by asymptotic scaling laws, using nonparametric multivariate functional estimation theory to extrapolate information from student models to teacher models.

3. We develop KCas for smoothing spline ANOVA models, design the associated algorithm, and establish consistency theory. Under the proposed subsampling scheme, KCas reduces the computational complexity of smoothing-parameter selection from \(O(n^3)\) to \(O(n^{\frac{3}{4}})\) while preserving statistical accuracy.

4. We illustrate the broader applicability of KCas beyond smoothing splines through kernel density estimation and deep learning hyperparameter transfer.

5. Through extensive simulations and real-data experiments, we show that KCas achieves substantial computational savings and can even outperform full-sample tuning in some settings.

2 Related Work↩︎

An idea related to KCas is self-distillation (SD), which is also developed to deploy complex models effectively. SD methods use the same architecture for both the teacher and student models, and facilitate the training procedure by letting the knowledge be transferred/exchanged among a group of models or within a single model [10], [17]–[22]. However, models in standard SD methods are still relatively large, and the model training procedure is accelerated and improved by distilling knowledge from itself. In this sense, KCas differs from SD by using information from some ‘actually small’ student models that are much easier to train. In scientific scenarios where the pilot study is needed to determine the design of extensive and detailed experiments, KCas can use the pilot data to construct student models to avoid wasting valuable data in the pilot study. Thus, KCas is highly desirable in these settings. Note that in [23], the authors also discussed the possibility of reversing the knowledge distillation procedure, but their methodology is still under the SD framework. [23] reverses the KD procedure as a motivating example for proposing the Teacher-free Knowledge Distillation (Tf-KD) framework. They prove the equivalence between KD and label smoothing regularization in a certain sense, and using this fact, Tf-KD lets a student model learn from itself or manually designed regularization distribution. Therefore, the student model in Tf-KD serves the purpose of regularization, while the student model in KCas serves the role of extracting information, and KCas amplifies this information to help the teacher model. Thus, our proposed KCas and the associated theories are significantly different from [23] and various self-distillation methods.

Further, [24] proposes a similar idea with KCas of distilling the knowledge from the student model to help train the teacher model. They first train a noisy student model on the available labeled data, which is then used to generate pseudo-labels for the unlabeled data. Subsequently, the teacher model is then trained on the combined set of labeled and pseudo-labeled data. On the contrary, our research focuses on addressing situations where the model training process requires an exceptionally high computational burden or when completing the necessary computations becomes impractical. Therefore, [24] addresses scenarios where the training model has sufficient computational resources but insufficient data, whereas our work focuses on situations where there is abundant data but limited computational resources available for training the model.

From a methodological perspective, KCas can be viewed as a strategy for efficient hyperparameter selection. Recent advancements in hyperparameter selection techniques continue to refine the efficacy of model selection. General approaches of hyperparameter selection include the classic cross-validation [25], [26], criterion-based methods [27], [28], and a comprehensive review of recent developments can be found in [29]. In the context of smoothing spline models, a series of methods based on generalized cross-validation (GCV) has been developed [30]–[34]. In this paper, we propose to conduct hyperparameter selection via the idea of knowledge cascade.

3 Preliminaries↩︎

In this section, we provide the necessary background of nonparametric functional estimation in reproducing kernel Hilbert space before introducing our proposed method. To estimate a function of interest \(\eta\) on a generic domain \(\mathcal{X}\), we consider the nonparametric penalized loss functional, \[\label{PLik} PL = L(\eta)+ {\lambda} J(\eta),\tag{1}\] where \(L(\eta)\) is the goodness-of-fit (loss) functional and \(J(\eta)\) is the smoothness (penalty) functional. The smoothing parameter \(\lambda\) controls the trade-off between the smoothness of \(\eta(x)\) and its fidelity to the data.

The estimation of function \(\eta\) on the product domain \(\mathcal{X}=\prod_{j=1}^{d} \mathcal{X}_{j}\) has long been a crucial problem in statistical learning, with numerous proposed methods over the years [35]–[39]. However, most of them have been challenged by the curse of dimensionality, as the estimation of multivariate functions is intrinsically difficult. To overcome this challenge, one effective approach is to decompose multivariate functions using techniques similar to the classical analysis of variance (ANOVA) decomposition and its associated notions of the main effect and interaction [35], [40]–[43]. In this functional ANOVA model, higher-order interactions are often excluded in practical estimation to control model complexity; excluding all interactions yields the popular additive models. On the product domain \(\mathcal{X}=\prod_{j=1}^{d} \mathcal{X}_{j}\), the function \(\eta\) can be decomposed as a sum of a constant term, one-dimensional functions (main effects), two-dimensional functions (two-way interactions), and so on, as in the following decomposition, \[\label{eq:decomposition} \eta(x)=\eta\left(x_{\langle 1\rangle}, \ldots, x_{\langle d\rangle}\right)=\eta_{\emptyset}+\sum_{j} \eta_{j}\left(x_{\langle j\rangle}\right)+\sum_{j<k} \eta_{j, k}\left(x_{\langle j\rangle}, x_{\langle k\rangle}\right)+\ldots,\tag{2}\] with the constant in \(\eta_{\emptyset}\), the main effects in \(\eta_{j}\), the two-way interactions in \(\eta_{j, k}\), etc. Higher-order interactions are eliminated to ease the curse of dimensionality.

By adding the smoothness penalty \(J(\eta)\) to \(L(\eta)\) in loss (1 ), we consider the space \(\mathcal{H} \subseteq\{\eta: J(\eta)<\) \(\infty\}\) in which \(J(\eta)\) is a square semi-norm with a finite-dimensional null space \(\mathcal{N}_{J}=\{\eta: J(\eta)=0\}\). To assist analysis and computation, a metric and geometry should be defined in this space, and the loss (1 ) needs to be continuous in \(\eta\) under this metric. Since the reproducing kernel Hilbert spaces (RKHSs) are well suited for this purpose, we consider the space \(\mathcal{H}\) as an RKHS with the continuous evaluation \([x] f=f(x)\), reproducing kernel \(R(\cdot, \cdot)\), a non-negative definite function satisfying \(\langle R(x, \cdot), f(\cdot)\rangle=f(x), \forall f \in \mathcal{H}\), where \(\langle\cdot, \cdot\rangle\) is the inner product in \(\mathcal{H}.\) The existence of the minimizer in RKHS is guaranteed by [44]. Details are discussed in Appendix 8.

Density Estimation. Consider the situation that we have independently identically distributed (i.i.d.) data points \(x_i,\;i=1,\cdots n\), from an underlying data distribution \(p(x)\) on a bounded domain \(\mathcal{X}=\prod_{j=1}^{d} \mathcal{X}_{j}\). We aim to estimate \(p(x)\) based on observations \(x_i\). For the nonparametric setting, a naive maximum likelihood density estimation is meaningless without any nonintrinsic constraint, since it will fit a sum of delta function spikes at the sample points \(x_i\), which is obviously not an appealing estimate when the domain \(X\) is continuous. Thus, a penalized likelihood estimate (PLE) is a good candidate. Two intrinsic constraints come from the definition of a probability density that \(p(\cdot)\ge 0\) and \(\int_{\mathcal{X}} p dx=1\). Since these two constraints are difficult to handle directly in computation, a common approach [45], [46] is to estimate the log-density \(\eta(\cdot)\), which is free of the constraints through the transformation \[\label{eq:den95trans} p(x)= \frac{e^{\eta(x)}}{\int_{\mathcal{X}} e^{\eta(x)} d x}.\tag{3}\] [46] proposed and studied the theoretical properties of the penalized likelihood over a Hilbert space \(\mathcal{H}\): \[\label{PLE} -\frac{1}{n} \sum_{i=1}^{n} \eta\left(x_{i}\right) + \log \int_{\mathcal{X}} e^{\eta (x)} d x+ \frac{\lambda}{2} J(\eta).\tag{4}\] Consider the exponential family with the densities of the form \[\label{eq:exp} f(y \mid x)=\exp \{\frac{y \eta(x)-h(\eta(x))}{a(\phi)}+c(y, \phi)\},\tag{5}\] where \(a(\cdot)>0\), \(h\), and \(c\) are known functions, \(\eta(x)\) is the regression function via the link \(\eta\), and \(\phi\) is the parameter that is independent of \(x\). Observing \(Y_i \mid x_i \sim f(y\mid x_i)\), \(i=1,\cdots ,n\), we estimate \(\eta(x)\) via the penalized likelihood functional \[\label{eq:PL95reg} -\frac{1}{n} \sum_{i=1}^{n}\left\{Y_{i} \eta\left(x_{i}\right)-h\left(\eta\left(x_{i}\right)\right)\right\}+\frac{\lambda}{2} J(\eta),\tag{6}\] where the term \(c(y,\phi)\) is dropped as it is independent of \(\eta(x)\), and \(a(\phi)\) is absorbed into \(\lambda\).

4 Methodology↩︎

Nonparametric penalized estimation of the function of interest \(\eta\) is a general question in lots of fields [47], [48]. However, the computational burden of training complex models limits the applicability of many existing methods to large datasets. To tackle this challenge, we propose a reverse version of knowledge distillation, named knowledge cascade, by cascading the knowledge learned from a student model (trained on a small sample) to the teacher model (trained on a large sample). Specifically, we illustrate our method in the context of nonparametric functional estimation, with two important cases: density estimation and regression functional estimation.

4.1 Minimizer of the Penalized Loss Functional↩︎

We first introduce the computation for the minimizer of the general penalized loss functional (1 ). Consider a tensor sum decomposition of the RKHS \({\cal H}= {\cal N}_J \oplus {\cal H}_{J}\). Let \(\{\phi_v\}_{v=1}^M\) be a basis of \(\mathcal{N}_J = \{\eta: J(\eta) = 0\}\) and \(R_J\) be the reproducing kernel in \(H_J\). Taking the ANOVA decomposition (2 ) into consideration, the RKHS \({\cal H}_J\) can be further decomposed into \({\cal H}_J= \oplus_{\beta=1}^g {\cal H}_{\beta}\) with the reproducing kernel \(R_J=\sum_{\beta=1}^g \theta_\beta R_\beta\), where \(R_{\beta}\) is the reproducing kernel in \({\cal H}_{\beta}\). Here the \(\theta_{\beta}\)’s are an extra set of smoothing parameters adjusting the contribution of the corresponding components. The minimizer of loss (1 ) is achieved in the tensor product reproducing kernel Hilbert space \(\mathcal{H}\) with the smoothness penalty \(J(\eta)=J(\eta, \eta)=\sum_{\beta=1}^g \theta_\beta^{-1}(\eta, \eta)_\beta\), where \((\eta, \eta)_\beta\) are inner products in \(\mathcal{H}_{\beta}\) with reproducing kernels \(R_\beta\). Without loss of generality, we define \(J(\eta)\) with tensor-product cubic splines throughout the paper. To ease the notation, we give an example of a tensor product cubic spline on \([0,1]^2\) in Appendix C, and please refer to [16] for the explicit forms of \(\{R_\beta\}_{\beta=1}^g\) and \(J(\eta)\).

According to the Kimeldorf-Wahba representer theorem [44], [49], [50], the minimizer of loss (1 ) has the following form \[\label{eq:eta95minimizer} \eta(x) = \sum_{v=1}^M d_v\phi_v(x)+\sum_{i=1}^{n} c_i R_J(x_i,x) = \boldsymbol{\varphi}(x)^{T} \boldsymbol{d}+\boldsymbol{\xi}(x)^{T} \boldsymbol{c},\tag{7}\] where \(\boldsymbol{d}= (d_1,\cdots,d_M)^T\), \(\boldsymbol{c}= (c_1,\cdots,c_n)^T\) are unknown coefficients, \(\boldsymbol{\varphi}(x)=(\phi_1, \cdots , \phi_M)^T\), \(\boldsymbol{\xi}(x)=(R_J(x_i,\cdot), \cdots , R_J(x_n,\cdot))^T\) are vectors of functions. Taking advantage of the representer theorem, the infinite-dimensional optimization problem is transferred into a finite-dimensional one, thereby facilitating the estimation.

For the density estimation problem, plugging the representer of \(\eta(x)\) (7 ) into the penalized likelihood of density estimation (4 ), the estimation reduces to the minimization of \[\label{eq:loss95den} -\frac{1}{n} \mathbf{1}^{T}(Q \mathbf{c}+S \mathbf{d})+\log \int \exp \left(\boldsymbol{\varphi}(x)^{T} \mathbf{d}+\boldsymbol{\xi}(x)^{T} \mathbf{c}\right) \mathrm{d} x+ \frac{\lambda}{2} \mathbf{c}^{T} Q \mathbf{c},\tag{8}\] where \(Q\) is \(n\times n\) with the \((i,j)\)th entry \(R_J(x_i, x_j)\) and \(S\) is \(n\times M\) with the \((i,v)\)th entry \(\phi_v(x_i)\). Similarly, the minimization of the penalized likelihood functional (6 ) for regression can be achieved via the minimizer of: \[\label{eq:loss95reg} \frac{1}{n} \left( \tilde{\boldsymbol{Y}} - S \mathbf{d} -Q \mathbf{c} \right)^T W \left( \tilde{\boldsymbol{Y}} - S \mathbf{d} -Q \mathbf{c} \right) + \frac{\lambda}{2} \mathbf{c}^{T} Q \mathbf{c},\tag{9}\] where \(W\) is the weight matrix, and the explicit form of \(\tilde{\boldsymbol{Y}}\) and \(W\) can be found in Appendix 9. Fixing smoothing parameters \(\lambda\) and \(\theta\), we can estimate the coefficients \(\boldsymbol{d}\) and \(\boldsymbol{c}\) in (8 ) or (9 ) using Cholesky decomposition [51] or Newton-Raphson method [16], [45].

To make the estimation practical, a critical aspect is selecting appropriate values for \(\lambda\) and \(\theta\) that result in satisfactory performance, since the solution of the penalized loss functional (1 ) is sensitive to \(\lambda\) and \(\theta\) [16], [35]. Smoothing parameters control the trade-off between the smoothness of \(\eta(x)\) and its fidelity to the data. Selecting a large smoothing parameter will lead to oversmoothing, while a small one will lead to undersmoothing.

One of the most efficient criteria for selecting the smoothing parameters is generalized cross-validation (GCV) [34], [52], which achieves the selection via cross-validation targeting the Kullback–Leibler (KL) loss. GCV consists of two main steps: (i) minimizing the KL loss with respect to \(\lambda\) for fixed \(\theta\); (ii) updating \(\theta\) according to the updated \(\lambda\). However, the computational cost of tuning the parameters, particularly \(\lambda\), can be prohibitively high in high-dimensional settings. With all \(S\) smoothing parameters tunable, the above iterative algorithm takes \(O(Sn^3)\) flops per iteration and needs tens of iterations to converge [34]. Here the number of smoothing parameters \(S\) increases as the number of multi-way interaction terms grows. In particular, \(\eta(x)=\eta\left(x_{\langle 1\rangle }, \ldots, x_{\langle d\rangle }\right)=\eta_{\emptyset}+\sum_{j} \eta_{j}\left(x_{\langle j\rangle }\right)+\sum_{j<k} \eta_{j, k}\left(x_{\langle j\rangle }, x_{\langle k\rangle }\right)\), the ANOVA decomposition model (2 ) truncated at two-way interactions contains \(S=d+\frac{3}{2}d(d-1)\) smoothing parameters. Considering the computational burden, in the case of a particularly large sample size, it is impractical to apply GCV on the full sample, i.e., train the teacher model directly, to accomplish the regression and density estimation tasks using the smoothing spline ANOVA model. To tackle this problem, we propose the knowledge cascade (KCas) method, which enables the determination of the smoothing parameters without incurring a significant computational burden during estimation.

4.2 Knowledge Cascade↩︎

In KCas, we aim to let student models learn the smoothing parameters through optimization, with a significantly lower computational burden compared to the teacher, and then transfer the smoothing parameters to teacher models. We first illustrate the definitions of the student and teacher models in our context of nonparametric functional estimation. The teacher model is a complex model trained on the full sample, and the student model is a simple model trained on a subset of the sample with a size of \(b\). In (7 ), we can see that the number of kernels, i.e., \(R_J(x_i,x)\), equals the sample size \(n\). Since the number of kernels highly affects the representation power of the model, the model complexity differs significantly for large and small sample sizes and thus distinguishes the student and teacher models. We first illustrate the simplest version of KCas in the general regression model with additive noise, \[\label{eq:ref} Y_i = \eta(x_i) + \epsilon(x_i),\tag{10}\] where \(\epsilon(\cdot)\) is the white noise process satisfying \(\text{E}(\epsilon(x_i))=0\), \(\text{E}(\epsilon(x_i)\epsilon(x_j))=\sigma^2\) if \(x_i=x_j\), \(\text{E}(\epsilon(x_i)\epsilon(x_j))=0\) otherwise. Define Hilbert space \({\cal H}^{(m)}\) by \[\begin{align} {\cal H}^{(m)}=\left\{\eta: \eta^{(v)}\right. &\text{absolutely continuous for } v=0,1, \ldots, m-1, \eta^{(m)} \in \mathcal{L}_{2}[0,1],\\ &\left.\eta^{(v)}(0)-\eta^{(v)}(1)=0 \text{ for } v=0,1, \ldots, m-1 \right\} , \end{align}\] where \(\eta^{(v)}= \frac{\mathrm{d}^v \eta}{ \mathrm{d} x^v}\), \(\mathcal{L}_{2}[0,1]=\left\{f: \int_0^1 f^2 d x<\infty\right\}\) is the Hilbert space formed by all square-integrable functions, and \(m\) is a constant indicating the order of smoothness of \({\cal H}^{(m)}\). For smoothing splines in \({\cal H}^{(m)}\), the optimal smoothing parameter \(\lambda\), ignoring \(o(1)\) terms, is \[\label{eq:lambda} \lambda = C n^{-\frac{2 m}{2 mp+1}},\tag{11}\] where \(C\) is an unknown constant depending on unknown function \(\eta\) [53] and \(p\in [1,2]\) indicates different additional smoothness conditions [33], [53], [54]. The estimation of \(C\) is infeasible since it depends on the unknown true function \(\eta\). However, KCas can infer the information of \(C\) from a well-trained subsample model (student) and apply it to the full data model (teacher). Specifically, notice that the asymptotically optimal \(\lambda_{\mathrm{GCV}}^{\mathrm{sub}}\) when sample size equals \(b\) is \[\lambda_{\mathrm{GCV}}^{\mathrm{sub}}(b) = Cb^{-\frac{2 m}{2 mp+1}},\] for the same \(C\). We estimate the optimal \(\lambda_{\mathrm{GCV }}^{\mathrm{sub}}\) on the subsample to infer the constant \(C\), and then employ the same \(C\) for the full data [55]. That is, we have the following estimation of \(\lambda\) by KCas, \[\label{KCas} \lambda_{\mathrm{KCas}}^{\mathrm{full}}(n;b) = \lambda_{\mathrm{GCV }}^{\mathrm{sub}}(b) (\frac{n}{b})^{-\frac{2m}{2mp+1}}.\tag{12}\] Since the smoothing parameters are used to determine the proportion of the smoothness penalty on different terms in (2 ) and this proportion should be stable over different sample sizes, we directly use the optimal \(\theta_{\mathrm{GCV}}^{\mathrm{sub}}(b)\) in the full sample. We then generalize the estimator (12 ) from the regression model with additive noise (10 ) to a wide range of penalized likelihood estimation problems, including the nonparametric regression in the exponential family and density estimation.

We propose to use KCas to transfer the knowledge of smoothing parameters in (4 ) or (6 ) to the teacher model. The KCas algorithm is summarized in the following Algorithm 1. In the first step, we apply uniform sampling to get a subsample \(X_b\), and our experiments show that uniform sampling can achieve good performance. Other more dedicated sampling methods can also be applied to improve the performance further [56]–[58]. The total number of operations required for each iteration is generally \(\frac{4}{3}n^3+O(n^2)\), in which the selection of the smoothing parameter takes the major burden. The GCV algorithm for the student model takes \(O(Sb^3)\) flops per iteration, and thus KCas algorithm reduces the computational cost from \(O(Sn^3)\) to \(O(Sb^3)\). According to the justifications in [42], [59]–[61], it is sufficient to take subsample size \(O(n^{\frac{2}{9}})\) to maintain the performance of smoothing spline estimation, and we thus take a slightly larger subsample size \(b = O(n^{\frac{1}{4}})\) for practical use in our algorithm. Consequently, our proposed KCas method, as detailed in Algorithm 1, substantially lowers the burden of the smoothing parameter estimation process to \(O(Sn^{\frac{3}{4}})\). This underscores the pivotal role of KCas.

In practice, hyperparameters \(m\) and \(p\) need to be selected properly. For the univariate setting, the commonly used smooth level is \(m=2\), which means the penalty is \(J(\eta, \eta)=\int_{0}^{1}\left(\eta^{(2)}\right)^{2} \mathrm{~d} x\) on \([0,1]\). For the multivariate setting, we suggest applying the commonly used tensor product cubic spline, which has \(2-\epsilon < m < 2, \forall \epsilon >0\) [44]. When \(\eta^{(2)}\) is square-integrable, we have \(p=1\), and when \(\eta^{(4)}\) is square-integrable, we have \(p=2\). The selection of \(p\) can vary across datasets. In this work, we take \(m=2\) and \(p=2\) empirically.

We should mention that GCV is typically used for Gaussian-type regression. For regression with responses from exponential families which results in a nonquadratic loss (6 ), the computation of \(\eta_\lambda\) requires a bunch of iterations even with fixed smoothing parameters and thus results in a more expensive computational cost. To address this issue, we use the well-accepted generalized approximate cross-validation (GACV) method [31], [62] or its variants to reduce the computational burden. Similarly, a variant of GCV [40] will be used to effectively select the smoothing parameter. To make the notation concise, variants of GCV, suitable for different problems, are collectively referred to as ‘GCV’.

4.3 Theoretical Analysis↩︎

In this section, we present the theoretical properties of the smoothing parameters \(\lambda\) selected according to Algorithm 1. For notational simplicity, in the following, we will use \(\lambda\) to represent all the smoothing parameters, not just the \(\lambda\) in front of \(J(\eta)\). All proofs for this section are relegated to Appendix 12.

Denote by \(\widehat{\eta}\) the estimate through the minimization of (4 ) and \(\eta_0\) the true function to be estimated, the asymptotic convergence rate based on the selected smoothing parameter \(\lambda_{\mathrm{KCas}}^{\mathrm{full}}(n;b)\) is established through Theorem 1.

With Theorem 1, the consistency is ensured, and the convergence rate is specified for the estimation \(\widehat{\eta}\) in density estimation based on KCas. Denote by \(\widehat{\eta}\) the estimate through the minimization of (6 ), the following Theorem 2 further demonstrates the theoretical properties of KCas in the context of nonparametric regression.

Remark. Note that the rates (?? ) and (?? ) are free of dimension \(d\), which is because we adopt the dimensionless approach as in Chapter 9 in [16], which essentially incorporates assumptions on the smoothness to get rid of the dimension in the rates. Specifically, Condition [cond:eigenvalue] implicitly puts an assumption on the smoothness of \(\eta\), which makes the error rates (?? ) and (?? ) to be free of dimension \(d\). The condition becomes more stringent as dimension \(d\) becomes larger. The results can be adapted to the case with dimension using arguments in [63].

4.4 Knowledge Cascade in Kernel Density Estimation↩︎

We next illustrate the KCas framework in the classical univariate kernel density estimation (KDE) problem. Let \(x_1,\ldots,x_n\) be i.i.d.observations from an unknown density \(f\) on \(\mathbb{R}\). For a symmetric kernel \(K\) and a bandwidth \(h>0\), the Gaussian kernel density estimator is \[\widehat{f}_h(x) = \frac{1}{n}\sum_{i=1}^n K_h(x-x_i) = \frac{1}{nh} \sum_{i=1}^n K\!\left(\frac{x-x_i}{h}\right),\] where we take \(K\) to be the standard normal density throughout [64]. For twice differentiable \(f\), the classical asymptotic mean integrated squared error (AMISE) theory yields an optimal bandwidth of the form \[h_{\mathrm{AMISE}} = \Bigg\{\frac{R(K)}{m_2(K)^2\,R\!\big(f''\big)}\Bigg\}^{1/5} n^{-1/5} = C\,n^{-1/5},\] where \(R(g)=\int g^2(x)\,dx\), \(m_2(K)=\int x^2 K(x)\,dx\), and \(C\) is an unknown constant depending on \(f\) and \(K\) [64], [65].

In KCas, we treat a KDE fitted on a small subsample as the student model and the full-sample KDE as the teacher model. The key idea is to estimate the constant \(C\) from a subsample of size \(b\) and then transfer it to the full data via the AMISE scaling. Specifically, the AMISE-optimal bandwidth for the subsample satisfies \[h_{\mathrm{AMISE}}^{\mathrm{sub}}(b) = C\,b^{-1/5},\] so that \(C = h_{\mathrm{AMISE}}^{\mathrm{sub}}(b)\,b^{1/5}\). Plugging this expression into the full-sample formula gives the KCas bandwidth \[\label{eq:KCas95KDE} h_{\mathrm{KCas}}^{\mathrm{full}}(n; b) = h_{\mathrm{AMISE}}^{\mathrm{sub}}(b)\!\left(\frac{n}{b}\right)^{-1/5}.\tag{13}\] In our implementation, \(h_{\mathrm{AMISE}}^{\mathrm{sub}}(b)\) is not computed by oracle access to \(f\) but by a data-driven selector applied to the subsample; the KCas bandwidth then extrapolates this choice to the full data. This is a direct KDE analogue of the smoothing-parameter scaling used in smoothing splines.

The same knowledge cascade principle extends naturally beyond KDE to a broad class of kernel-based methods. In many kernel procedures, such as kernel regression [66], local polynomial smoothing [67], and kernel-based covariance estimation [68], the optimal tuning parameters follow explicit sample-size dependent scaling laws derived from asymptotic theories. KCas exploits this structure by learning the problem-dependent constants from a small subsample and transferring them to the full sample through the corresponding scaling relations, thereby reducing computational cost while preserving statistical efficiency. This perspective suggests KCas as a general strategy for tuning-parameter transfer in kernel methods with known asymptotic rates.

4.5 Knowledge Cascade in Deep Learning↩︎

Although KCas is developed in the context of nonparametric multivariate functional estimation, the underlying principle extends naturally to settings in which computational constraints make full hyperparameter tuning difficult. In many deep learning applications, the cost of identifying stable and high-performing hyperparameter configurations, particularly the learning rate schedule, dominates the total training budget. Motivated by this challenge, we investigate a simple extension of KCas that uses a compact student network to guide the hyperparameter choices for a larger teacher network.

In this setting, the student model is trained on a reduced architecture and its optimal hyperparameters are obtained through standard tuning procedures. KCas then scales these hyperparameters to the teacher model using a batch-size dependent rule. This approach acknowledges empirical findings that learning rate schedules tuned on smaller networks often transfer to larger networks when batch sizes are adjusted appropriately. Consequently, KCas reduces hyperparameter search cost by allowing tuning to be performed once on the student model and then extrapolated to the teacher.

Mathematically, let \(\eta_{\text{ student }}\) denote the optimal learning rate for the student model, and let \(B_{\text{student }}\) and \(B_{\text{teacher }}\) be the respective batch sizes. The KCas learning rate for the teacher model is defined as \[\eta_{\text{ teacher }}=\eta_{\text{ student }} \cdot g\left(\frac{B_{\text{teacher }}}{B_{\text{student }}}\right) ,\] where \(g(\cdot)\) is a monotone scaling function. We consider two commonly used choices. The first is the linear rule \(g(r)=r\), inspired by the large-batch scaling observations in [69]. Since teacher architectures in our experiments are substantially more complex and often benefit from more conservative step sizes, we also examine the square-root rule \(g(r)=\sqrt{r}\). These two options illustrate how KCas can adapt scaling behavior to model complexity.

Although this extension does not rely on the RKHS-based theory that supports KCas in nonparametric functional estimation, it follows the same guiding philosophy: use a computationally efficient student model to extract stable hyperparameter information and propagate it upward to a more complex teacher model. As demonstrated in Section 6.2, this strategy provides a practical reduction in tuning cost while maintaining competitive predictive performance for modern deep learning architectures.

5 Simulation Study↩︎

We conduct simulation studies to evaluate the proposed KCas framework from two perspectives. First, we examine its performance in the main setting of this paper: smoothing spline ANOVA models for density estimation and nonparametric regression. These experiments assess whether smoothing parameters learned from a subsample-based student model can be effectively transferred to the full-sample teacher model while reducing computational cost. Second, we investigate KCas in KDE to illustrate the same student-to-teacher transfer principle in a different nonparametric setting where bandwidths follow classical asymptotic scaling laws.

5.1 Simulation 1: Density Estimation with Smoothing Splines↩︎

We first evaluate the proposed KCas method for smoothing splines on synthetic density estimation problems. We consider the following two data-generating scenarios.
Scenario 1: A \(d\)-dimensional Gaussian mixture model, \(\frac{1}{d} \sum_{i=1}^d \text{MVN}(e_i, I_d)\), where \(e_i\) is the vector with the \(i\)th entry being \(1\) and others being \(0\). We consider \(d=3,6\).
Scenario 2: A \(d\)-dimensional density is constructed by independently combining a \(5\)-dimensional Gaussian with mean zero and variance \(0.5(\mathbf{1}\mathbf{1}^T)+1.5I_d\), and the remaining \(d-5\) variables are i.i.d. from Unif\((0,1)\). We consider \(d=15,20\).

We evaluate the methods using the relative Kullback–Leibler (KL) divergence with respect to the benchmark method, defined as \[\operatorname{RelKL}\left(\widehat{p}, p^*, p\right) = \frac{D_{\mathrm{KL}}(\widehat{p}\Vert p)}{D_{\mathrm{KL}}(p^* \Vert p)}, \nonumber\] where \(p\) is the true distribution, \(\widehat{p}\) is the estimate from the method being evaluated, and \(p^*\) is the estimate from the benchmark method. Here \(D_{\mathrm{KL}}(q\Vert p)\) denotes the KL divergence of \(q\) relative to \(p\). A smaller relative KL divergence indicates better performance. In each replication, full-sample GCV is used as the benchmark, and the relative KL divergence is computed with respect to the corresponding GCV estimate from that replication. For visualization, we report log-transformed RelKL.

For the smoothing spline experiments, including Simulation 1 here and following Simulation 2, we compare KCas with the following baseline methods:

• GCV [40] on the full sample: Full-sample GCV uses all available observations for smoothing-parameter selection and is therefore used as the benchmark. All other methods are evaluated using relative performance measures with respect to full-sample GCV.

• GCV on subsample (SUB): GCV is performed on a randomly selected subsample to reduce computational cost.

• GCV in generalized additive models (GAM) [70]: GCV is performed on the full sample using a simplified additive model with only main effects.

• SKIP method [30]: SKIP accelerates computation by selecting an appropriate starting point and bypassing subsequent iterations. Since it often fails to converge in relatively high-dimensional density estimation settings, we include it only for the nonparametric regression problem.

• Order-based method (ORD) [32]: ORD directly sets the smoothing parameter to \(\lambda=n^{-\frac{2m}{2mp+1}}\), where \(n\) is the full sample size.

• Kernel density estimation (KDE): For density estimation, we further include the KDE method of [71] as a comparison for high-dimensional data.

For KCas, we use uniform sampling to select a subsample of size \(b=50n^{\frac{1}{4}}\). The full sample sizes are \(5,000\), \(10,000\), and \(20,000\). Following the discussion in Section 4.2, we set \(m=2\) and \(p=2\) in practice. All results are based on \(30\) replications.

Figure 2 illustrates the log-transformed RelKL of KCas and the competing methods. KCas performs favorably across the considered settings, with log-RelKL values close to or below zero in many cases. Negative values indicate that KCas outperforms the full-sample GCV benchmark in that replication. This suggests that smoothing parameters learned from the student model can effectively guide the teacher model and, in some settings, may also reduce the effect of overly adaptive full-sample tuning.

5.2 Simulation 2: Nonparametric Regression with Smoothing Splines↩︎

For the nonparametric regression problem, we consider the model \[y_i \sim \text{Ber}\left(\frac{\exp\left(\eta\left(x_i\right)\right)}{1+\exp\left(\eta\left(x_i\right)\right)}\right), \nonumber\] where \(\text{Ber}(p)\) denotes the Bernoulli distribution with probability \(p\), \(x_{i}=\left(x_{i\langle 1\rangle}, \ldots, x_{i\langle d\rangle}\right)^{\mathrm{T}}\) is the \(d\)-dimensional predictor for the \(i\)th observation, with each entry independently drawn from Unif\((0,1)\). \(\eta\) is the nonparametric function determining the success probability in the Bernoulli trial, and \(y_{i} \in \{0,1\}\) is the response variable for the \(i\)th observation.

We evaluate the methods by the relative MSE, defined by \[\text{RelMSE}\left(\widehat{\eta}, \eta^*, \eta\right) = \frac{\sum_{i=1}^{n}\left\{\widehat{\eta}\left(x_{i}\right)-\eta\left(x_{i}\right)\right\}^{2}}{\sum_{i=1}^{n}\left\{\eta^{*}\left(x_{i}\right)-\eta\left(x_{i}\right)\right\}^{2}} , \nonumber\] where \(\eta\) is the true function, \(\widehat{\eta}\) is the estimation by the method being evaluated, and \(\eta^*\) is the benchmark method. We compare KCas with the same baseline methods as in Simulation 1.

We consider two different scenarios with different dimensions. We report log-transformed RelMSE for plotting clearness, and each RelMSE is computed based on a full GCV baseline for that particular replication.
Scenario 1: \[\begin{align} \eta_{m 1}(x) = \sum_{i=1}^3 g_1\left(x_{\langle i \rangle}\right) + g_2\left(x_{\langle 1 \rangle}, x_{\langle 2 \rangle}\right) + g_2\left(x_{\langle 1 \rangle}, x_{\langle 3 \rangle}\right) + g_3\left(x_{\langle 1 \rangle}, x_{\langle 2 \rangle}, x_{\langle 3 \rangle}\right). \end{align}\]

Scenario 2: \[\begin{align} \eta_{m 2}(x) &= \sum_{i=1}^3 i g_1\left(x_{\langle i \rangle}\right) + \sum_{i=4}^6 i g_{5}\left(x_{\langle i \rangle}\right) + \sum_{i=7}^9 g_{4}\left(x_{\langle i \rangle}\right) + \\ & \sum_{i=1}^3 \sum_{j>i}^4 3i g_2\left(x_{\langle i \rangle}, x_{\langle j \rangle}\right) + 6 g_2\left(x_{\langle 5 \rangle}, x_{\langle 6 \rangle}\right) + 8 g_6\left(x_{\langle 7 \rangle}, x_{\langle 8 \rangle}\right) + 10 g_3\left(x_{\langle 1 \rangle}, x_{\langle 2 \rangle}, x_{\langle 3 \rangle}\right). \end{align}\]

The explicit forms of the functions \(g_i\) are provided in Appendix 13.1. Scenario 1 is a well-established setting for nonparametric multivariate functional estimation in RKHS [34], [35], [41], [55]. We consider two cases, with \(d=3\) and \(d=6\) predictors. When \(d=3\), all predictors contribute to \(\eta_{m1}\) and hence to the response \(y\). When \(d=6\), the last three predictors are irrelevant to \(\eta_{m1}\), representing a setting with redundant covariates. Scenario 2 is a more complex high-dimensional setting, where we consider \(d=15\) and \(d=20\).

Figure 3 displays the log-transformed RelMSE relative to full-sample GCV. In Scenario 1 with \(d=3\), KCas and SKIP show comparable performance and outperform the other methods. Their median log-RelMSE values are close to, and sometimes below, zero, indicating performance comparable to or better than full-sample GCV. When \(d=6\), KCas remains close to full-sample GCV, whereas SKIP performs substantially worse, with median log-RelMSE values greater than \(1\). This is likely because SKIP uses the starting values introduced by [30] as the final estimate; when \(\eta\) becomes more complex, such starting values may be far from the optimum, leading to inaccurate estimation. In Scenario 2, similar patterns are observed for both \(d=15\) and \(d=20\), with KCas achieving the best overall performance.

To empirically compare computational efficiency, we report the running time under the most challenging simulation setting, Scenario 2 with \(d=20\), in Table 1. KCas substantially reduces computation time compared with full-sample GCV while maintaining comparable estimation accuracy. It is important to note that while methods such as GAM, SUB, and SKIP require less computation time than KCas, they do not achieve the same level of accuracy.

5.3 Simulation 3: Kernel Density Estimation↩︎

We next conduct an additional simulation to evaluate KCas in classical kernel density estimation. This experiment is separate from the smoothing spline simulations above and uses bandwidth-selection methods as baselines. The goal is to examine whether information learned from a small subsample can be transferred to the full-sample KDE estimator through the classical \(n^{-1/5}\) bandwidth scaling rule.

We consider six univariate benchmark densities widely used in KDE studies [64], [65]: (1) standard normal density \(\mathcal{N}(0,1)\); (2) symmetric bimodal Gaussian mixture \(0.6\,\mathcal{N}(-1,1^2)+0.4\,\mathcal{N}(2,0.5^2)\); (3) lognormal density \(\text{Lognormal}(0,1)\); (4) trimodal Gaussian mixture \(0.4\,\mathcal{N}(-2,0.5^2)+0.3\,\mathcal{N}(0,0.6^2)+0.3\,\mathcal{N}(3,0.8^2)\); (5) gamma density \(\text{Gamma}(2,1.5)\); and (6) skewed bimodal Gaussian mixture \(0.7\,\mathcal{N}(-1,1^2)+0.3\,\mathcal{N}(3,0.5^2)\). For each density, we generate samples of size \(n \in \{200, 500, 1000, 2000, 5000\}\). For a given density and \(n\), we draw \(50\) independent replications and, for each replication, compute several bandwidth selectors and their associated KDEs.

We compare the following methods:

• Oracle AMISE: the theoretical AMISE-optimal bandwidth \(h_{\mathrm{AMISE}}\) obtained by numerically evaluating \(R\!\big(f''\big)\) from the known \(f\). This serves as a lower bound on the achievable mean integrated squared error (MISE).

• Improved Sheather-Jones (ISJ): the diffusion-based plugin selector of [72], which refines the original Sheather–Jones solve-the-equation method [73] and is known to perform well across a wide range of densities.

• Least-squares cross-validation (LSCV): the classical least-squares cross-validation bandwidth [74], [75], which minimizes an unbiased estimate of the integrated squared error.

• Silverman’s rule of thumb: the normal-reference bandwidth \(h_{\mathrm{S}} = 0.9\min\{\widehat{\sigma},\mathrm{IQR}/1.34\}\) \(n^{-1/5}\) [64], where \(\widehat{\sigma}\) is the sample standard deviation and \(\mathrm{IQR}\) is the sample interquartile range, defined as the difference between the third and first sample quartiles, \(\mathrm{IQR} = Q_{3} - Q_{1}\).

• KCas-Silverman: a KCas estimator based on the Silverman rule computed on a uniform subsample of size \(b=200\), i.e.\(h_{\mathrm{KCas},S}(n;b)\) obtained by plugging \(h_{\mathrm{AMISE}}^{\mathrm{sub}}(b) = h_{\mathrm{S}}^{\mathrm{sub}}(b)\) into 13 .

• KCas-ISJ: a KCas estimator based on the ISJ selector computed on a subsample of size \(b=200\), i.e.\(h_{\mathrm{KCas},\mathrm{ISJ}}(n;b)\) from 13 with \(h_{\mathrm{AMISE}}^{\mathrm{sub}}(b)\) taken as the subsample ISJ bandwidth.

• KCas-CV: a KCas estimator based on least-squares cross-validation (LSCV) computed on a subsample of size \(b=200\), i.e.\(h_{\mathrm{KCas},\mathrm{CV}}(n;b)\) obtained by plugging the subsample LSCV bandwidth \(h_{\mathrm{CV}}^{\mathrm{sub}}(b)\) into 13 .

In each method, we fit a Gaussian KDE on an equally spaced grid of \(4096\) points covering the effective support of each density. The integrated squared error (ISE) for a replication is approximated via the trapezoidal rule, \[\mathrm{ISE}(\widehat{f}_h) \approx \sum_{k} \big\{\widehat{f}_h(x_k) - f(x_k)\big\}^2 \,\Delta x,\] and we summarize performance using the median ISE and interquartile range across the \(50\) replications. To facilitate comparison, we also report the ratio of the median ISE to that of the oracle AMISE bandwidth for each density and \(n\), as well as log–log slopes of median bandwidth versus sample size to check the \(n^{-1/5}\) scaling.

Figure 4 displays the distribution of ISE values across replications for all selectors, stratified by density and sample size. As expected, the oracle AMISE bandwidth achieves the smallest median ISE in every setting and provides a lower bound for the other procedures. Among implementable full-sample selectors, the diffusion-based ISJ bandwidth is consistently the strongest competitor, with median ISEs very close to the oracle across all \(n\) and for both simple (normal, gamma) and more complex (bimodal, trimodal, skewed) densities, in line with previous empirical studies [65], [72], [73]. Silverman’s rule-of-thumb performs competitively for the unimodal, approximately Gaussian densities, but becomes noticeably suboptimal for heavier-tailed and highly multimodal shapes, where its normal-reference assumption leads to oversmoothing. The LSCV selector exhibits the largest variability, especially for smaller sample sizes, with wider ISE distributions and occasional outlying bandwidths, reflecting its well-known finite-sample instability [65], [74], [75].

The KCas variants largely inherit the strengths and weaknesses of their corresponding base selectors while operating on only a small subsample. For the normal and gamma densities, KCas-Silverman attains median ISE ratios on the order of \(1.05\)–\(1.20\) relative to the oracle across the range of \(n\), indicating only a modest loss in efficiency despite learning bandwidths from \(b=200\) observations. For multimodal and skewed densities, KCas-ISJ is more adaptive and typically dominates KCas-Silverman, with median ISE curves that nearly overlap those of the full-sample ISJ selector. Across all densities, the variability of KCas-ISJ remains close to that of ISJ itself, suggesting that the additional randomness induced by subsampling does not substantially degrade performance. KCas-CV performs reasonably well in smooth unimodal settings, particularly for the normal density, but is generally less competitive for heavier-tailed and multimodal densities and exhibits larger variability, reflecting the weaker finite-sample stability of cross-validation based selectors. Overall, KCas-Silverman and KCas-ISJ consistently outperform KCas-CV, benefiting from more structured statistical foundations and stronger theoretical guarantees.

From the computational perspective, the cost of KCas is dominated by fitting the student model on \(b=200\) points, rendering the bandwidth selection cost essentially independent of \(n\) once \(b\) is fixed. Taken together, these results demonstrate that, in the KDE setting, KCas provides a practical mechanism for transferring the adaptivity of sophisticated selectors such as ISJ to large-sample regimes, while maintaining optimal scaling and substantially reducing the computational burden of hyperparameter tuning.

Simulation results across various settings and methods show that KCas can effectively transfer knowledge of smoothing parameters from the student model on the subsample to the teacher model on the full sample, consistent with our theoretical analysis.

6 Real Data Analysis↩︎

We apply KCas to real datasets in two settings: smoothing-spline-based nonparametric estimation and deep learning hyperparameter transfer. Specifically, we consider density estimation and nonparametric regression on benchmark datasets, as well as image classification on CIFAR-10.

6.1 Nonparametric Estimation↩︎

We first evaluate KCas on smoothing-spline-based nonparametric estimation tasks. The density estimation study uses four benchmark datasets, and the nonparametric regression study uses five benchmark datasets. Details of the datasets are provided in Appendix 14. Each dataset is randomly split into 80% training and 20% testing sets. Features are scaled using a min–max transformation fitted on the training set and then applied to the testing set. All relative metrics in this subsection are computed with respect to full-sample GCV.

Density estimation. We compare KCas with GAM, SUB, ORD, and KDE, in density estimation. Since there is no ground truth for the density function, we evaluate the performance based on the average log-likelihood on the test set, as suggested by previous studies [76], [77]. We consider the ANOVA decomposition of \(\eta\) including all main effects and all two-way interactions. The model terms are selected using the model diagnosis suggested by [78]. Our results, presented in Table 2, demonstrate that KCas outperforms all four benchmark methods in terms of log-likelihood across all data sets. Notably, on the ESC and MFCC datasets, KCas even outperforms GCV on the full sample, both in terms of log-likelihood and computational time.

Nonparametric regression. We compare KCas with GAM, SUB, ORD, and SKIP, in nonparametric regression for exponential families. Since we do not know the underlying probability of each data point, as suggested by [56], we calculate RelMSE over GCV. The main effect and interaction terms are selected using the model diagnosis suggested by [78]. Table 3 shows that KCas outperforms all four benchmark methods in terms of MSE. Although KCas is not the fastest among the methods, it is faster than the full sample estimator in all studies, while obtaining the best performance among the comparison methods.

6.2 Image Classification with Deep Learning↩︎

We next apply KCas to deep learning hyperparameter transfer on the CIFAR-10 dataset [79]. CIFAR-10 contains \(50,000\) training images and \(10,000\) test images from ten evenly represented object classes: airplane, automobile, bird, cat, deer, dog, frog, horse, ship, and truck.

Model Structure. We consider two student–teacher pairs. MobileNetV2 with 2.2M parameters [80] and ResNet18 with 11.7M parameters [81] serve as student models, while ResNet50 with 25.6M parameters serves as the teacher model. All architectures are adapted to the CIFAR-10 resolution by removing max-pooling and replacing the initial layers with a \(3\times3\) convolution. Standard dataset normalization is applied.

KCas Hyperparameter Search. To obtain student configurations for transfer, each student model is tuned by grid search over \(12\) hyperparameter combinations: batch sizes \({128,256,512}\), learning rates \({0.05,0.1}\), and weight-decay values \({5\mathrm{e}{-4},1\mathrm{e}{-3}}\). Students are trained for \(200\) epochs with SGD using cosine decay and a 5% warmup. The best configuration is selected using a 10% holdout validation set. KCas then scales the selected learning rate to the teacher model.

Baselines. We compare KCas with three baselines. (1) Student model: the accuracy of the tuned student model. (2) Cookbook: hand-selected hyperparameters commonly recommended for the teacher architecture. (3) Retune: a full grid search for the teacher over \(36\) combinations: batch sizes \({256,512,1024}\), learning rates \({0.02,0.05,0.1,0.2}\), and weight-decay values \({1\mathrm{e}{-4},5\mathrm{e}{-4},1\mathrm{e}{-3}}\). This baseline represents the standard but computationally expensive approach to teacher hyperparameter tuning.

Results. For each method, we record test accuracy and training time, averaged over three repetitions.

Figure 5 shows that KCas can effectively transfer learning-rate information from smaller networks to a larger teacher network, achieving teacher performance competitive with full hyperparameter re-tuning. Panel (a) compares KCas with the baselines. Both KCas variants, using linear and square-root scaling, substantially improve over the student models and reach accuracy levels close to those of fully re-tuned teachers for both MobileNetV2 and ResNet18. The cookbook rule provides moderate gains but consistently falls short of the proposed KCas approaches. An additional pattern in panel (a) is that teacher performance is higher when ResNet18 is used as the student rather than MobileNetV2. This likely reflects the stronger standalone performance of ResNet18, which allows it to provide more informative hyperparameters for transfer. The results suggest that the student model should not be too weak in practice. Although MobileNetV2-based transfer still outperforms the cookbook baseline, a stronger student yields clearer benefits for the teacher.

Panel (b) directly compares the two KCas scaling rules across all batch-size choices, rather than only the best setting. Linear and square-root scaling achieve similar accuracy overall, while the square-root rule shows slightly greater stability across repetitions. Together with panel (a), these results suggest that square-root scaling tends to perform better in this experiment.

Table 5 reports the computational time. Compared with full re-tuning of the teacher model, which requires more than \(120\) hours, KCas uses only a small fraction of the computational budget. The savings are more pronounced with a smaller student model, since the total cost of KCas depends on the student-training time. In practical scenarios where a trained student model is already available, the additional cost of KCas reduces to training a single teacher model, making it comparable to the cookbook baseline. Although KCas is slightly more expensive than the cookbook approach when the student must be trained from scratch, it consistently achieves better performance. These results illustrate that KCas offers a favorable balance between efficiency and predictive performance for deep learning hyperparameter transfer.

7 Discussion↩︎

In this paper, we propose Knowledge Cascade (KCas), a student-to-teacher knowledge transfer framework developed primarily for nonparametric functional estimation in RKHS. In this setting, a student model fitted on a small subsample guides smoothing-parameter selection for the full-sample teacher model, reducing tuning cost while retaining strong statistical performance. Our simulation and real-data results show that KCas compares favorably with several computationally efficient alternatives and can even outperform full-sample GCV in some settings. Beyond the main RKHS setting, we further illustrate the same student-to-teacher transfer principle through kernel density estimation and deep learning hyperparameter transfer. These examples suggest that KCas is not limited to smoothing spline models, but can provide a broader strategy for scalable model development whenever useful information learned from a smaller model can be transferred to a larger one. More broadly, KCas offers a new perspective on knowledge distillation by showing that small models can sometimes guide, rather than merely approximate, larger models. It remains an open question whether such knowledge cascades can be constructed without the support of asymptotic or empirical scaling rules.

When and why to use KCas. KCas is most useful when a smaller student model can learn information that remains informative for developing a larger teacher model. The key requirement is not that the student approximates the teacher perfectly, but that it captures transferable information that can reduce the cost of training or tuning the teacher. Such information may include smoothing parameters, bandwidths, training configurations, or other structural features that are expensive to learn directly at full scale. When this information follows a suitable asymptotic or empirical scaling relationship, KCas provides a way to extract it from a low-cost student model and propagate it to the teacher.

KCas is not intended to replace the teacher with a smaller approximation. Instead, the student serves as a computationally efficient guide for building the teacher. In this paper, we instantiate this idea mainly through hyperparameter transfer. For example, in smoothing spline ANOVA models, KCas learns smoothing parameters from a subsample-based student model and transfers them to the full-sample teacher model through an asymptotic scaling rule. Thus, KCas avoids repeated full-scale tuning while still fitting the final model on all observations. More broadly, the same principle suggests a general strategy for reducing the cost of teacher-model development whenever useful information can be learned cheaply from a smaller model and reliably transferred to a larger one.

KCas can outperform the full sample estimator. In multiple scenarios in our simulation and real data analysis, we have observed an interesting phenomenon that KCas, based on a random subsample, can outperform the full sample GCV estimator. One potential explanation is that the statistical theory underlying KCas helps make hyperparameter transfer more efficient. Although this finding appears to conflict with the traditional intuition, recent studies [82]–[86], have observed similar phenomena that models trained on subsamples sometimes can achieve better performance in the context of deep models. Further, a related paper [87] observed such a phenomenon for simpler settings under empirical risk minimization, and they aim to develop some theoretical justification for it. Although they work under different model settings, it is an interesting future direction to explore this phenomenon theoretically and empirically in nonparametric estimation or deep learning settings.

Connection to broader statistical methodology. The KCas principle may also be relevant to other statistical methods whose effective complexity changes with sample size. Examples include basis-expansion methods, wavelet estimators, series regression, and other multiresolution procedures. In these settings, quantities such as the number of basis terms, resolution levels, thresholding parameters, or regularization strengths often need to be adjusted as the sample size increases. This makes them natural candidates for KCas-type extensions, since a student model fitted at a smaller scale may help estimate problem-specific quantities that guide the corresponding full-scale estimator. Establishing suitable scaling rules and theoretical guarantees for these broader classes of methods remains an important direction for future work.

Code availability. The source code is available at: https://github.com/LuyangFang/KCas.

8 Existence of the Minimizer↩︎

The following theorem guarantees the existence of the minimizer of (1 ) in RKHS.

When \(L(\eta)\) is the negative log-likelihood function, it is usually convex in \(\eta\). The quadratic functional \(J(\eta)\) is also convex [45]. A minimizer of \(L(\eta)\) is unique in \(\mathcal{N}_{J}\) if the convexity is strict in it, which is often the case. Thus, the solution for Equation (1 ) exists in most cases.

9 Minimizer of the Penalized Loss Functional↩︎

Consider exponential family distributions with densities of the form \[f(y \mid x)=\exp \{\frac{y \vartheta(x)-b(\vartheta(x))}{a(\phi)}+c(y, \phi)\},\] where \(a>0, b\), and \(c\) are known functions, \(\vartheta(x)\) is the canonical parameter dependent on a covariate \(x\), and \(\phi\) is either known or considered as a nuisance parameter that is independent of \(x\). Fixing the smoothing parameters, the penalized likelihood functional (6 ) is strictly convex in \(\eta\). Thus, given the current \(\tilde{\eta}\), the Newton iteration can be used to update \(\tilde{\eta}\) by the minimizer of the penalized weighted least square functional

\[\frac{1}{n} \left( \tilde{\boldsymbol{Y}} - S \mathbf{d} -Q \mathbf{c} \right)^T W \left( \tilde{\boldsymbol{Y}} - S \mathbf{d} -Q \mathbf{c} \right) + \frac{\lambda}{2} \mathbf{c}^{T} Q \mathbf{c},\] where \(\tilde{Y}_i=\tilde{\eta}\left(x_i\right)- \frac{\tilde{u}_i}{\tilde{w}_i}\), \(\tilde{w}_i=\ddot{b}\left(\tilde{\eta}\left(x_i\right)\right), \text{and } \tilde{u}_i=-Y_i+\dot{b}\left(\tilde{\eta}\left(x_i\right)\right)\). Here \(\tilde{w}_i\) is the \(i\)th diagonal element of the matrix \(W\), \(\tilde{Y}_i\) is the \(i\)th element of \(\tilde{\boldsymbol{Y}}\), \(\dot{b}\) and \(\ddot{b}\) are the first and second derivatives of the function \(b\).

10 Example of the Tensor Product Cubic Spline on \([0,1]^2\)↩︎

For the univariate \(\eta\) on \(\mathcal{X}\), the most popular choice of the smoothness penalty \(J(\eta)\) is \[J(\eta,\eta) = \int_{\mathcal{X}} \left(\eta^{(m)}\right)^2 \mathrm{~d} x,\] where \(\eta^{(m)}=\mathrm{d}^m \eta / \mathrm{d} x^m\). A cubic estimator of the minimizer of (1 ) is obtained by setting \(m=2\). This idea can be extended to multivariate settings. Consider the ANOVA decomposition (2 ). \(\mathcal{H}\) can be decomposed into the space of constants, the spaces of main effects, and the corresponding spaces of interaction terms lying in the tensor product space of the interacting main-effect spaces. For the two-dimensional problem, one has the following space decomposition in each variable ([16], section 2.3): \[\begin{align} \left\{\eta: \eta^{(2)} \in \mathcal{L}_2[0,1]\right\}= & \{\eta: \eta \propto 1\} \oplus\left\{\eta: \eta \propto k_1\right\} \\ & \oplus\left\{\eta: \int_0^1 \eta \mathrm{~d} x=\int_0^1 \eta^{(1)} \mathrm{d} x=0, \eta^{(2)} \in \mathcal{L}_2[0,1]\right\} \\ = & \mathcal{H}_{00} \oplus \mathcal{H}_{01} \oplus \mathcal{H}_1, \end{align}\] where \(k_1(x)=x-0.5\). The space of constant terms is \(\mathcal{H}_{00\langle1\rangle} \otimes \mathcal{H}_{00\langle2\rangle}\); the space of main effects is spanned by \(\mathcal{H}_{00\langle1\rangle} \otimes\left(\mathcal{H}_{01\langle2\rangle} \oplus \mathcal{H}_{1\langle2\rangle}\right)\) and \(\mathcal{H}_{00\langle 2\rangle} \otimes\left(\mathcal{H}_{01\langle 1\rangle} \oplus \mathcal{H}_{1\langle1\rangle}\right)\); and the subspace \(\left(\mathcal{H}_{01\langle 1\rangle} \oplus \mathcal{H}_{1\langle 1\rangle}\right) \otimes\) \(\left(\mathcal{H}_{01\langle2\rangle} \oplus \mathcal{H}_{1\langle2\rangle}\right)\) spans the space of interactions. Let \(\mathcal{H}_{\nu, \mu}=\mathcal{H}_{\nu\langle 1\rangle} \otimes \mathcal{H}_{\mu\langle 2\rangle}\) for \(\nu, \mu=00,01,1\), with inner products \((\eta, \eta)_{\nu, \mu}\) and reproducing kernels \(R_{\nu, \mu}=R_{\nu\langle 1\rangle} R_{\mu\langle 2\rangle}\), using the tensor product cubic spline, we have \[\begin{align} J(\eta, \eta)= & \theta_{1,00}^{-1}(\eta, \eta)_{1,00}+\theta_{00,1}^{-1}(\eta, \eta)_{00,1} \\ & +\theta_{1,01}^{-1}(\eta, \eta)_{1,01}+\theta_{01,1}^{-1}(\eta, \eta)_{01,1}+\theta_{1,1}^{-1}(\eta, \eta)_{1,1} , \end{align}\] and the null space of \(J(\eta, \eta)\) is \[\mathcal{N}_J=\mathcal{H}_{00,00} \oplus \mathcal{H}_{01,00} \oplus \mathcal{H}_{00,01} \oplus \mathcal{H}_{01,01} .\]

11 Regularity Conditions↩︎

We define the quadratic functional representing the mean square error of the estimator \(\widehat{\eta}\) in estimating the target function \(\eta_0\) on the domain \(\mathcal{X}\) as \[V\left(\widehat{\eta}-\eta_0 \right)=\int_{\mathcal{X}}\left\{\widehat{\eta}-\eta_0(x)\right\}^{2} f(x) \mathrm{d} x,\] where \(f(x)\) is the marginal density of \(x\). We now state some regularity conditions for Theorem 1 and Theorem 2.

When condition [cond:continuous] is satisfied, that is, \(V\) is completely continuous with respect to \(J\) and hence to \(V+J\), there exist eigenvalues \(\lambda_{\nu}\) and the corresponding eigenfunctions \(\psi_{\nu}\) such that \[\begin{align} &V\left(\psi_{\nu}, \psi_{\mu}\right)=\lambda_{\nu} \delta_{\nu, \mu}, \text{ and } \\ &(V+J)\left(\psi_{\nu}, \psi_{\mu}\right)=\delta_{\nu, \mu}, \end{align}\] where \(\delta_{\nu, \mu}\) is the Kronecker delta and \(1 \geq \lambda_{\nu} \downarrow 0\); see [88], [46].

Write \(\phi_{\nu}=\lambda_{\nu}^{-\frac{1}{2}} \psi_{\nu}\). It follows that \[\begin{align} V\left(\phi_{\nu}, \phi_{\mu}\right)&=\delta_{\nu, \mu},\\ J\left(\phi_{\nu}, \phi_{\mu}\right)&=\rho_{\nu} \delta_{\nu, \mu}, \end{align}\] where \(0 \leq \rho_{\nu}=\lambda_{\nu}^{-1}-1\). We refer to \(\rho_{\nu}\) as the eigenvalues of \(J\) with respect to \(V\) and to \(\phi_{\nu}\) as the associated eigenfunctions. A Fourier series expansion of \(\eta_0\) satisfying \(J(\eta_0)<\infty\) is \(\eta_0=\sum_{\nu} \eta_{\nu,0} \phi_{\nu}\), where \(\eta_{\nu,0}=V\left(\eta_0, \phi_{\nu}\right)\) are the Fourier coefficients.

Condition [cond:bound] requires a uniform bound for the fourth moments of \(\phi_{\nu}(X)\).

Condition [cond:equi] requires the equivalence of the information in \(B_0\).

12 Proofs of Main Results↩︎

Note that it has been proved rigorously that the optimal smoothing parameter \(\lambda(b)\) has the form \(C b^{-\frac{r}{rp+1}}\) under some exponential regression problems, such as regression with Gaussian-type responses and periodic splines [33], [53], [54]. In such cases, with the fact that \(rp>1\), as \(b \rightarrow \infty\), \[\begin{align} \lambda(b)&= C b^{-\frac{r}{rp+1}} \rightarrow 0, \text{ and} \\ b \lambda(b)^{\frac{2}{r}} &= b C^{\frac{2}{r}} b^{-\frac{2}{rp+1}} \\ &= C^{\frac{2}{r}}b^{\frac{rp-1}{rp+1}} \rightarrow \infty, \end{align}\] that is, \(\lambda(b) \rightarrow 0\) and \(b (\lambda(b))^{\frac{2}{r}} \rightarrow \infty\) is naturally satisfied. In some cases, such as the density estimation problems and more general exponential-family settings, we impose the stated assumptions on \(\lambda_{\mathrm{GCV}}^{\mathrm{sub}}(b)\); the numerical results support their validity.

We replace \(\lambda(b)\) with \(\lambda_{\mathrm{GCV}}^{\mathrm{sub}}(b)\) chosen by GCV since it is infeasible to determine \(\lambda(b)\) with the unknown function \(\eta_0\). Theoretical results [54], [89] have shown that \(\lambda_{\mathrm{GCV}}^{\mathrm{sub}}(b)\) is a good estimator of \(\lambda(b)\), with \(\frac{L(\lambda_{\mathrm{GCV}}^{\mathrm{sub}}(b))}{L(\lambda(b))}=1+o_p(1)\). Thus, it is natural to extend the assumption \(\lambda_{\mathrm{GCV}}^{\mathrm{sub}}(b) \rightarrow 0\) and \(b (\lambda_{\mathrm{GCV}}^{\mathrm{sub}}(b))^{\frac{2}{r}} \rightarrow \infty\) to the general regression problems with responses from exponential families. The numerical results also support this assumption.

13 Simulation Details↩︎

13.1 Nonparametric regression↩︎

Scenario 1: Let \[\begin{align} \eta_{m 1}(x) = \sum_{i=1}^3 g_1\left(x_{\langle i \rangle}\right) + g_2\left(x_{\langle 1 \rangle}, x_{\langle 2 \rangle}\right) + g_2\left(x_{\langle 1 \rangle}, x_{\langle 3 \rangle}\right) + g_3\left(x_{\langle 1 \rangle}, x_{\langle 2 \rangle}, x_{\langle 3 \rangle}\right), \end{align}\]

Scenario 2: Let \[\begin{align} \eta_{m 2}(x) &= \sum_{i=1}^3 i g_1\left(x_{\langle i \rangle}\right) + \sum_{i=4}^6 i g_{5}\left(x_{\langle i \rangle}\right) + \sum_{i=7}^9 g_{4}\left(x_{\langle i \rangle}\right) + \\ & \sum_{i=1}^3 \sum_{j>i}^4 3i g_2\left(x_{\langle i \rangle}, x_{\langle j \rangle}\right) + 6 g_2\left(x_{\langle 5 \rangle}, x_{\langle 6 \rangle}\right) + 8 g_6\left(x_{\langle 7 \rangle}, x_{\langle 8 \rangle}\right) + 10 g_3\left(x_{\langle 1 \rangle}, x_{\langle 2 \rangle}, x_{\langle 3 \rangle}\right), \end{align}\] where
\(g_1(x)=10^6 x^{11}(1-x)^6\);
\(g_2\left(x_{\langle 1 \rangle}, x_{\langle 2 \rangle}\right)=\exp \left(3 x_{\langle 1 \rangle} x_{\langle 2 \rangle}\right)\);
\(g_3\left(x_{\langle 1 \rangle}, x_{\langle 2 \rangle}, x_{\langle 3 \rangle}\right)=\) \(\frac{15 \sin \left(2 \pi x_{\langle 1 \rangle}\right)}{2-\sin \left(2 \pi x_{\langle 2 \rangle} x_{\langle 3 \rangle}\right)}\);
\(g_4(x)=10^4 x^{3}(1-x)^{10}\);
\(g_5(x)=15x sin(15x)\);
\(g6(x) = \frac{a p_1}{\pi \sigma_{1} \sigma_{2}} \exp \left\{-\frac{\left(x_{\langle 1\rangle}-0.2\right)^{2}}{\sigma_{1}^{2}}-\frac{\left(x_{\langle 2\rangle}-0.3\right)^{2}}{\sigma_{2}^{2}}\right\} +\frac{a p_2}{\pi \sigma_{1} \sigma_{2}} \exp \left\{-\frac{\left(x_{\langle 1\rangle}-0.7\right)^{2}}{\sigma_{1}^{2}}-\frac{\left(x_{\langle 2\rangle}-0.8\right)^{2}}{\sigma_{2}^{2}}\right\}-b,\) with \(\sigma_{1}=0.3\), \(\sigma_{2}=0.4\), \(p_1=0.625\), \(p_2=0.375\), and \(a=b=4.2\).

14 Datasets↩︎

14.1 Datasets for Density Estimation↩︎

• CD14: Transcriptions in CD14 single cells. The data contains the abundance information of 13 proteins in 2,096 cells. The data set is available from [90].

• AReM: Activity Recognition system based on Multisensor data fusion Data Set. The dimension is 6 and the sample size is 42,240. The time-domain features including 3 mean values and 3 standard deviations were collected from the multisensor system during a period of time. The data set is available at UCI Machine Learning Repository [91].

• ESC: Embryonic Stem Cell from Mouse [92]. The data concerns mouse embryonic stem cell gene expression and transcription factor association strength. The 4 features that describe the scores of TFAS with KLF4, NANOG, OCT4, and SOX2 of 1,027 genes are used for density estimation. The data set is available at CRAN in gss package.

• MFCC: Anuran Calls (MFCCs) Data Set. The data is extracted from syllables of anuran (frogs) calls, including 22 variables with a sample size of 7,195. The data set is available at UCI Machine Learning Repository [91].

All the continuous variables in these datasets are scaled through a min-max normalization.

14.2 Datasets for Nonparametric Regression↩︎

• SUSY: Supersymmetric Dataset [93]. The dataset contains one response and 18 kinematic features \(x_{\langle 1\rangle},\dots,x_{\langle 18\rangle}\). The full sample size is 5,000,000, and about \(54.24\%\) of the responses in the data are from the background process. We consider the full sample GCV as the gold standard, but it is not affordable to compute GCV on the full sample. Therefore, we randomly select a subsample of size \(20,000\) and consider this sample to be the “full sample” to compute the GCV, and only use this part of data to conduct all the following analysis in the paper. The data set is available at UCI Machine Learning Repository [91].

• WFRN: Wall-Following Robot Navigation Dataset [94]. The data is a robot navigating through the room following the wall using 24 ultrasound sensors with a sample size of 19,735. The data set is available at UCI Machine Learning Repository [91].

• OCUP: Room Occupancy Estimation Dataset [95]. The experimental testbed for occupancy estimation was deployed in a room. The setup consisted of 7 sensor nodes and one edge node in a star configuration with the sensor nodes transmitting data to the edge every 30s using wireless transceivers. The dataset contains one response and 15 features. The full sample size is 10,129. The data set is available at UCI Machine Learning Repository [91].

• SHILL: Shill Bidding Dataset [96]. This is a dataset with a large number of eBay auctions of a popular product. The dataset contains one response and 12 features. The full sample size is 6,321. The data set is available at UCI Machine Learning Repository [91].

• CIFAR-10: The CIFAR-10 dataset [79] consists of a training set of \(50,000\) examples and a test set of \(10,000\) examples. Each example in the dataset is a \(32\times32\) color image, spanning 10 different classes of objects such as animals and vehicles. These classes include airplanes, cars, birds, cats, deer, dogs, frogs, horses, ships, and trucks, each equally represented in the dataset.

  For this dataset, we use a pre-trained convolutional neural network to extract the features from the raw images first. This neural network model consists of two main components: a convolutional layer block and a fully-connected layer block. The convolutional block comprises two sets of convolutional layers with batch normalization and max pooling layers for feature extraction from input images. The fully-connected block contains four linear layers with ReLU activations. The third linear layer with 20 nodes serves as a feature extraction layer, providing a compressed representation of the input features for downstream tasks.

  For this experiment, we convert CIFAR-10 into a binary logistic regression task by setting \(Y=1\) for the car class and \(Y=0\) for all other classes. The MSE is computed on the test set between the binary responses and the fitted conditional probabilities, and the relative MSE is normalized by the MSE of the full-sample estimator.

All the continuous predictors in these datasets are scaled through a min-max normalization. For all the datasets, we apply the cosine diagnostics [16] first for the identifiability and the practical significance of the fitted terms, in order to avoid overfitting and overinterpreting.

References↩︎

---

1. Equal contribution.↩︎

2. Corresponding authors.↩︎