Oversampling Divide-and-conquer for Response-skewed Kernel Ridge Regression

We consider the problem of calculating large-scale kernel ridge regression (KRR) estimates in a nonparametric regression model. Although the theoretical properties of the KRR estimator are well-understood [8, 50, 38], in practice, the computation of KRR estimates may suffer from a large computational burden. In particular, for a sample of size $N$, it requires $O(N^{3})$ computational time to calculate a KRR estimate using the standard approach, as will be detailed in Section 2. Such a computational cost is prohibitive when the sample size $N$ is considerable. The divide-and-conquer approach has been implemented pervasively to alleviate such computational burden [51, 52, 44, 45, 46]. Such an approach randomly partitions the full sample into $k$ subsamples of equal sizes, then calculates a local estimate on an independent local processor (also called a local node) for each subsample. The local estimates are then averaged to obtain the global estimate. The divide-and-conquer approach reduces the computational cost of calculating KRR estimates from $O(N^{3})$ to $O(N^{3}/k^{2})$. Such savings may be substantial as $k$ grows.

Despite algorithmic benefits, the success of the divide-and-conquer approach highly depends on the assumption that the subsamples can well represent the observed full sample. Nevertheless, this assumption cannot be guaranteed in many real-world applications, where the response variable $Y$ may have a highly skewed distribution. Specifically, the random variable $Y$ has a highly skewed distribution $f_{Y}$ if $f_{Y}$ is nearly zero inside a large region $A_{Y}\subset\mathbb{R}$. Problems of this type arise in high energy physics, Bayesian inference, financial research, biomedical research, environmental data, among others [28, 1, 12]. In these applications, $Y\in A_{Y}$ are the responses that occur with reduced frequency. However, such responses are often of more interest as they tend to have a more widespread impact. For example, $Y\in A_{Y}$ may represent a rare signal for seismic activity or stock market fluctuations. Overlooking such signals could resulting in a substantial negative impact on society either economically or in terms of human casualties.

Recall that in the classical divide-and-conquer approach, a random subsample is processed on every local node. Under the aforementioned rare-event scenarios, such a subsample could fail to have observations selected from the potentially informative region $A_{Y}$. The local estimate based on such a subsample thus is very likely to overlook the region $A_{Y}$. Averaging these local estimators could lead to unreliable estimations and predictions over these informative regions. A synthetic example in Fig. 1 illustrates the scenario that the response is highly skewed. In this example, the one-dimensional sample (gray points) is uniformly generated from $[0,1]$. The response variable $Y$ has a heavy-tailed distribution, as illustrated by the histogram. The classical dacKRR method is used to estimate the true function (gray curve). We observe that almost all the local estimates (blue curves) Averaging these local estimates thus results in a global estimate (red curve) with unacceptable estimation performance. Such an observation is due to the fact that the subsample in each node overlooks the informative region over the distribution of $Y$ (the peaks) and thus fails to provide a successful global estimate.

Our contributions. To combat the obstacles, we develop a novel response-adaptive partition approach with oversampling to obtain large-scale kernel ridge regression estimates, especially when the response is highly skewed. Although the oversampling technique is widely used for addressing discrete label imbalance, extending such a technique to the continuous dacKRR setting is nontrivial. We bridge the gap by providing both theoretical and practical guidance on how to effectively over-sample the observations. Different from the classical divide-and-conquer approach, such that each observation is allocated to only one local node, we propose to allocate some informative data points to multiple nodes. Theoretically, we show the proposed estimates have a smaller risk than the classical dacKRR estimates under mild conditions. Furthermore, we show the number of strata $l$ of the response vector and the number of nodes $k$ regulate the trade-off between the computational time and the risk of the proposed estimator. In particular, when $k$ is well-selected, a larger $l$ is associated with a longer computational time as well as a smaller risk. Such results are novel to the best of our knowledge. Our theoretical findings are supported by both simulated and real-data analyses. In addition, we show the proposed method is not specific to the dacKRR method and has the potential to improve the estimation of other nonparametric regression estimators under the divide-and-conquer setting.

Model setup. Let ${\mathcal{H}}$ be a reproducing kernel Hilbert space (RKHS). Such a RKHS is induced by the reproducing kernel $K(\cdot,\cdot):\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$, which is a symmetric nonnegative definite function, and an inner product $\langle\cdot,\cdot\rangle_{\mathcal{H}}$ satisfying $\langle g(\cdot),K(\bm{x},\cdot)\rangle_{\mathcal{H}}=g(\bm{x})$ for any $g\in{\mathcal{H}}$. The RKHS ${\mathcal{H}}$ is equipped with the norm $||g||_{\mathcal{H}}=(\langle g(\cdot),g(\cdot)\rangle)^{1/2}$. The well-known Mercer’s theorem states that, under some regularity conditions, the kernel function $K(\cdot,\cdot)$ can be written as $K(\bm{x},\bm{z})=\sum_{j=1}^{\infty}\mu_{j}\phi_{j}(\bm{x})\phi_{j}(\bm{z})$, where $\{\mu_{j}\}_{j=1}^{\infty}\geq 0$ is a sequence of decreasing eigenvalues and $\{\phi_{j}(\cdot)\}_{j=1}^{\infty}$ is a family of orthonormal basis functions. The smoothness of a function $g\in{\mathcal{H}}$ is characterized by the decaying rate of the eigenvalues $\{\mu_{j}\}_{j=1}^{\infty}$. The major types of such decaying rates include the finite rank type, the exponentially decaying type, and the polynomially decaying type. The representative kernel function of these three major types are the polynomial kernel $K(\bm{x},\bm{z})=(1+\bm{x}^{\intercal}\bm{z})^{r}$ ($r$ is an integer), the Gaussian kernel $K(\bm{x},\bm{z})=\mbox{exp}(-||\bm{x}-\bm{z}||^{2}/\sigma^{2})$ ($\sigma>0$ is the scale parameter and $||\cdot||$ is the Euclidean norm), and the kernel $K(x,z)=1+\mbox{min}(x,z)$, respectively. We refer to [14, 36] for more details.

Consider the nonparametric regression model

where $y_{i}\in\mathbb{R}$ is the response, $\bm{x}_{i}\in\mathbb{R}^{d}$ is the predictors, $\eta$ is the unknown function to be estimated, and $\{\epsilon_{i}\}_{i=1}^{N}$ are i.i.d. normal random errors with zero mean and unknown variance $\sigma^{2}<\infty$. The kernel ridge regression estimator aims to find a projection of $\eta$ into the RKHS ${\mathcal{H}}$. Such an estimator $\widehat{\eta}$ can be written as

Here, the regularization parameter $\lambda$ controls the trade-off between the goodness of fit of $\widehat{\eta}$ and the smoothness of it. A penalized least squares framework analogous to Equation (2) has been extensively studied in the literature of regression splines and smoothing splines [40, 13, 23, 16, 10, 33, 48, 37, 47, 24, 49, 29].

The well-known representer theorem [40] states that the minimizer of Equation (2) in the RKHS ${\mathcal{H}}$ takes the form $\widehat{\eta}(\bm{x})=\sum_{i=1}^{N}\beta_{i}K(\bm{x}_{i},\bm{x}).$ Let $\bm{y}=(y_{1},\ldots,y_{N})^{\intercal}$ be the response vector, $\bm{\beta}=(\beta_{1},\ldots,\beta_{N})^{\intercal}$ be the coefficient vector and ${\mathbf{K}}$ be the kernel matrix such that the $(i,j)$-th element equals $K(\bm{x}_{i},\bm{x}_{j})$. These coefficient vector $\bm{\beta}$ can be estimated through solving the minimization problem as follows,

It is known that the solution of such a minimization problem has a closed form

where the regularization parameter $\lambda$ can be selected based on the cross-validation technique, the Mallows’s Cp method [27], or the generalized cross-validation (GCV) criterion [41].

Although the solution of the minimization problem (3) has a closed-form, the computational cost for calculating the solution using Equation (4) is of the order $O(N^{3})$, which is prohibitive when the sample size $N$ is considerable. In this paper, we focus on the divide-and-conquer approach for alleviating such a computational burden. In recent decades, there also exist a large number of studies that aim to develop low-rank matrix approximation methods to accelerate the calculation of kernel ridge regression estimates [31, 42, 32, 26, 30]. These approaches are beyond the scope of this paper. In practice, one can combine the divide-and-conquer approach with the aforementioned methods to further accelerate the calculation.

Background of the divide-and-conquer approach. The classical divide-and-conquer approach is easy to describe. Rather than solving the minimization problem (3) using the full sample, the divide-and-conquer approach randomly partitions the full sample into $k$ subsamples of equal sizes. Each subsample is then allocated to an independent local node. Next, $k$ minimization problems are solved independently on each local node based on the corresponding subsamples. The final estimate is simply the average of all the local estimates. The divide-and-conquer approach has proven to be effective in linear models [4, 21], partially linear models [53], nonparametric regression models [51, 52, 19, 35, 11], principal component analysis [43, 6], matrix factorization [25], among others.

Algorithm 1 summarizes the classical divide-and-conquer method under the kernel ridge regression setting. Such an algorithm reduces the computational cost for the estimation from $O(N^{3})$ to $O(N^{3}/k^{2})$. The savings may be substantial as $k$ grows.

One difficulty in Algorithm 1 is how to choose the regularization parameter $\lambda_{j}$s. A natural way to determine the size of these regularization parameters is to utilize the standard approaches, e.g., the Mallows’s Cp criterion or the (GCV) criterion, based only on the local subsample $\SS_{j}$. It is known that such a simple strategy may lead to a global estimate that suffers from suboptimal performance [52, 45]. To overcome the challenges, there has been a large number of studies dedicated to developing more effective methods of selecting the local regularization parameters in the recent decade. For example, [52] proposed to select the regularization parameter $\lambda_{j}$ according to the order of the entire sample size $N$ instead of the subsample size $N/k$. Recently, [45] proposed the distributed GCV method to select a global optimal regularization parameter for nonparametric regression under the classical divide-and-conquer setting.

Motivation. To motivate the development of the proposed method, we first re-examine the oversampling strategy. Such a strategy is well-known in imbalanced data analysis, where the labels can be viewed as skewed discrete response variables. [17, 15, 18]. Classical statistical and machine learning algorithms assume that the number of observations in each class is roughly at the same scale. In many real-world cases, however, such an assumption may not hold, resulting in imbalanced data. Imbalanced data pose a difficulty for classical learning algorithms, as they will be biassed towards the majority group. Despite the rareness, the minority class is usually more important from the data mining perspective, as it may carry useful and important information. An effective classification algorithm thus should take such unbalances into account. To tackle the imbalanced data, the oversampling strategy supplements the training set with multiple copies of the minority classes and keeps all the observations in the majority classes to make the whole dataset suitable for a standard learning algorithm. Such strategy has proven to be effective for achieving more robust results [17, 15, 18].

Intuitively, the data with a continuous skewed response can be considered as imbalanced data. In particular, if we divide the range of the skewed response $\{y_{i}\}_{i=1}^{N}$ into $l$ equally-spaced disjoint intervals, denoted by ${\mathcal{Y}}_{1},\ldots,{\mathcal{Y}}_{l}$, then there will be a disproportionate ratio of observations in each interval. The intervals that contain more observations can be considered as the majority classes, and the ones that contain fewer observations can be considered as the minority classes. Recall that the classical divide-and-conquer approach utilizes a simple random subsample from the full sample to calculate the local estimate. Therefore, when the response variable is highly skewed, such a subsample could easily overlook the observations respecting the minority classes, resulting in unpleasant estimates.

Main algorithm. Inspired by the oversampling strategy, we propose to supplement the training data with multiple copies of the minority classes before applying the divide-and-conquer approach. In addition, analogous to the downsampling strategy, the subsample in each node should keep most of the observations from the minority classes. Let $|\SS|$ be the number of observations in the set $\SS$. Let $[\cdot]$ be the rounding function that rounds down the number to its nearest integer. The proposed method, called oversampling divide-and-conquer kernel ridge regression, is summarized in Algorithm 2.

We use a toy example to illustrate Algorithm 2 in Fig. 2. Suppose there are three nodes ($k=3$), and the histogram of the response is divided into three slices ($l=3$). The left panel of Fig. 2 illustrates the oversampling process. In particular, the observations are duplicated by certain times, such that the total number of observations respecting each slice is roughly equal to each other. The original observations and the duplicated observations are marked by tilted lines and horizontal lines, respectively. The right panel of Fig. 2 shows the next step. For each slice, the observations within such a slice are then randomly and evenly allocated to each of the nodes. The observations allocated to different nodes are marked by blue, red and purple, respectively. Finally, similar to the classical dacKRR approach, the local estimates are calculated and then are averaged to obtain the global estimate.

Implementation details and the computational cost. Analogous to Algorithm 1, a difficulty in Algorithm 2 is how to choose the regularization parameters $\lambda_{i}$’s. Throughout this paper, we opt to use the same approach to select $\lambda_{i}$’s as the one proposed in [52]. The number of slices $l$ can be set as a constant or be determined by Scott’s rule; see [34] for more details. Our empirical studies indicate that the performance of the proposed algorithm is robust to a wide range of choices of $l$. In addition, to avoid the data dependency within a node raised by oversampling, we conduct a “de-duplicating” process before calculating the local estimates.

Consider the total number of observations in the training set after the oversampling and de-duplicating processes. Let $\tilde{N}$ denote such a number. When the response variable has a roughly uniform distribution, few observations would be copied, and thus one has $\tilde{N}=O(N)$. In such a case, the computational cost of Algorithm 2 is at the same order as that of the classical divide-and-conquer method. In the most extreme case that the most majority class ${\mathcal{Y}}_{1}$ contains almost all the observations, it can be shown that $\tilde{N}\leq lN$. Therefore, the computational cost of Algorithm 2 is at the order of $O(l^{3}N^{3}/k^{2})$. Consequently, when $l$ is a constant, again, the computational cost of Algorithm 2 is at the same order of the classical divide-and-conquer method. Otherwise, when $l$ and $k$ go to infinity, for example, when $l$ is determined by Scott’s rule [34] ( $l=O(N^{1/3})$), and $k=O(l)$, the computational cost of Algorithm 2 becomes $O(N^{4}/k^{2})$. However, such a higher computational cost is associated with theoretical benefits, as will be detailed in the next section.

Main theorem. In this section, we obtain the upper bounds on the risk of the proposed estimator $\bar{\eta}$ and establish the main result in Theorem 4.1. We show our bounds contain a smaller asymptotic estimation variance term, compared to the bounds of the classical dacKRR estimator [52, 20]. The following assumptions are required to bound the terms in Theorem 4.1. Assumption 4 guarantees the spectral decomposition of the kernel via Mercer’s theorem. Assumption 4 is a regularity assumption on $\eta_{0}$. These two assumptions are fairly standard. Assumption 4 requires that the basis functions are uniformly bounded. [52] showed that one could easily verify this assumption for a given kernel $K$. For instance, this assumption holds for many classical kernel functions, e.g., the Gaussian kernel.

Let $\{\tilde{{\mathbf{X}}}_{i},\tilde{\bm{y}}_{i}\}$ be the subset with size of $\tilde{n}$ allocated to the $i$th node after de-duplicating, $i=1,\ldots,k$. We assume they are drawn i.i.d from an unknown probability measure $\mathcal{P}$. Let $\mathcal{P}_{X}$ denote marginal distribution of $\tilde{{\mathbf{X}}}_{i}$. For brevity, we further assume the sample size for each node is the same. The bounds on the risk of the proposed estimator $\bar{\eta}$ are given in Theorem 4.1, followed by the main corollary.

With the properly chosen regularization parameter, most of the commonly used kernels, e.g., the kernels with polynomial eigendecay, can satisfy this assumption for Corollary 5 [9]. More discussions can be found in the Supplementary Material. The squared bias term for the classical dacKRR estimator is at the order of $O(\lambda\left\lVert\eta_{0}\right\rVert^{2}_{\mathcal{H}})$, which is the same as the one in Equation (5) [52]. However, different from the variance term in Equation (5), which equals $O(\sigma^{2}d_{\lambda}/\tilde{N})$, the variance term for the classical dacKRR estimator is at the order of $O(\sigma^{2}d_{\lambda}/N)$. Our rate depends on the eigenvalues $\{\mu_{j}^{2}\}$ instead of the trace of kernel derived in [52]. This rate is tighter than that in [2] and is considered as the minimax rate.

Practical guidance on the selection of the parameter $l$ and $k$. We now compare the risk of the classical dacKRR estimator with the proposed one under three different scenarios. Furthermore, we provide some practical guidance on the selection of the number of slides $l$ and the number of nodes $k$. We denote the total sample size after oversampling as $\tilde{N}^{*}$.

First, consider the scenario that $N\approx\tilde{N}^{*}$ as shown in Fig. 3(a). In this case, oversampling is unnecessary, thus the proposed estimate $\bar{\eta}$ has the same risk as the classical dacKRR estimator.

Second, consider the scenario that the response variable has a slightly skewed distribution, i.e., one has $\tilde{N}^{*}=O(N)$ as shown in Fig. 3(b). Under this scenario, the total sample size after de-duplicating is also in the same order of the original sample size, i.e. $\tilde{N}=O(N)$. Corollary 5 indicates the proposed estimate $\bar{\eta}$ has the same risk as the classical dacKRR estimator.

Third, consider the scenario that the response variable has a highly skewed distribution. Under this scenario, the sample size in the “majorest” slide is in the same order of $N$, and the sample size in the “minorest” slide, denoted as $N_{l}$, is in the order of $o(N)$. In these cases, one has $\tilde{N}^{*}=O(lN)$ We then discuss under two cases, (1) when $N_{l}=o(N/l)$, especially, when $l$ is a constant; and (2) when $N_{l}=O(N/l)$. In case (1), the sample size $N$ after de-duplicating is in the same order of $N$. The risk of $\bar{\eta}$ then has the same order as the classical dacKRR estimator. In case (2), notice that $N_{l}=o(N)$, this indicates that $l$ should go to infinity as $N\to\infty$. In this scenario, we consider $k=O(N/N_{l})$ i.e., $k=O(l)$, then we have $\tilde{N}=O(\tilde{N}^{*})$. Thus equation (5) becomes

For example, when $l$ is determined by Scott’s rule, one has $l=O(N^{1/3})$, and we choose $k=O(l)$. Equation (6) indicates the proposed estimator has a smaller risk than the classical dacKRR estimator. Corollary 5 indicates that the parameters $l$ and $k$ regulate the trade-off between the computational time and the risk. Specifically, when $k$ is well-selected, a larger $l$ is associated with a longer computational time as well as a smaller risk.

To show the effectiveness of the proposed estimator, we compared it with the classical dacKRR estimator in terms of the mean squared error (MSE). We calculated the MSE for each of the estimators based on 100 replicates. In particular, $\mbox{MSE}=\sum_{i=1}^{100}\|\hat{\eta}^{(i)}-\eta_{0}\|^{2}/100)$, where $\eta_{0}$ is the true function and $\hat{\eta}^{(i)}$ represents the estimator in the $i$th replication, respectively. Through all the experiments in this section, we set the number of nodes $k=100$. Gaussian kernel function was used for the kernel ridge regression. We followed the procedure in [52] to select the bandwidth for the kernel and the regularization parameters. For the proposed method, we divided the range of $\{y_{i}\}_{i=1}^{N}$ into $l$ slices according to Scott’s rule [34].

Recall that in Algorithm 2, each observation is copied for certain times, such that the total number of observations in each slice is roughly equal to each other. A natural question is how does the number of such copies affects the performance of the proposed method. In other words, if a fewer number of copies are supplemented to the training data, would the proposed method perform better or worse? To answer this question, we induced a constant $\tau$ to control the oversampling size. Specifically, for the $j$th slice, $j=1,\ldots,l$, we duplicated the observations within such a slice for $([\tau|{\mathcal{Y}}_{1}|/|{\mathcal{Y}}_{j}|]-1)$ times, instead of $([|{\mathcal{Y}}_{1}|/|{\mathcal{Y}}_{j}|]-1)$ times in Algorithm 2. We set $\tau=\{1,1/2,1/3,1/4,1/5\}$. This procedure is equivalent to Algorithm 2 when $\tau=1$.

Let $\bm{1}^{d}$ be the $d$-dimensional vector, such that all elements are equal to $1$. We considered a function that has a highly skewed response, i.e.,

where $\bm{x}\in[0,1]^{d}$, and $\|\cdot\|$ represents the $\mathbb{L}_{2}$ norm. We then simulated the data from Model (1) with $N=\{2,4,6,8,10\}\times 10^{3}$, $d=1,2,4,6$, and two different regression function $\eta_{0}$’s,

Uni-peak: $\eta_{0}(\bm{x})=g(\bm{x},\bm{c})$, with $\bm{c}=0.4\times\bm{1}^{d}$;

Double-peak: $\eta_{0}(\bm{x})=g(\bm{x},\bm{c}_{1})+0.4g(\bm{x},\bm{c}_{2})$, with $\bm{c}_{1}=0.4\times\bm{1}^{d}$ and $\bm{c}_{2}=0.7\times\bm{1}^{d}$.

Figure 4 shows the MSE versus different sample size $N$ under various settings. Each column represents a different $d$, and each row represents a different $\eta_{0}$. The classical dacKRR estimator and the full sample KRR estimator are labeled as blue lines and black lines, respectively. The proposed estimators are labeled as red lines, and the proposed method with $\tau=1$, i.e., Algorithm 2, is labeled as solid red lines. The vertical bars represent the standard errors obtained from 100 replications. In Fig. 4, we first observed that the classical dacKRR estimator, as expected, does not perform well. We then observed that the proposed estimators perform consistently better than the classical estimator. Such an observation indicates that when the response is highly skewed, oversampling the observations respecting the minority values of the response helps improve the estimation accuracy. Finally, we observed that a larger value of $\tau$ is associated with a better performance. In particular, as the number of copies increases, the proposed estimator tends to have faster convergence rates. This observation supports Corollary 5, which states that a larger number of observations after the oversampling process is associated with a smaller estimation MSE.

Besides Scott’s rule ($l=O(N^{1/3})$), we also considered other rules, e.g., the Sturge’s formula ($l=O(\log_{2}N)$) [39] and the Freedman-Diaconis choice ($l=O(N^{1/2})$) [7]. Figure 5 compares the empirical performance of the proposed estimator w.r.t. different choices of $l$’s. For the setting of unskewed response (the first row), all methods share similar performance. For response-skewed settings (the second and the third row), all the three aforementioned rules yield better performance than fixed $l$. Among all, Scott’s rule gives the best results. The simulation results are consistent with Corollary 4.2, which indicates that the proposed estimator shows advantages over the classical one.

Besides the impact of $l$, we also studied the impact of $k$. In addition, we have also applied the proposed strategy to other nonparametric regression estimators, say smoothing splines, under the divide-and-conquer setting. Such simulation results are provided in the Supplementary Material. These results indicated that the performance of the proposed method is robust the choice of $k$ and has the potential to improve the estimation of other nonparametric regression estimators under the divide-and-conquer setting. It is also possible to consider unequally-spaced slicing methods in the proposed method; however, such extensions are beyond the scope of this paper.

We applied the proposed method to a real-world dataset called Melbourne-housing-market ¹¹1Data source https://www.kaggle.com/anthonypino/ melbourne-housing-market.. The dataset includes housing clearance data in Melbourne from the year 2016 to 2017. Each observation represents a record for a sold house, including the total price and several predictors. Of interest is to predict the price-per-square-meter (range from 20 to 30000) for each sold house using the longitude and the latitude of the house, and the distance from the house to the central business district (CBD). Such a goal can be achieved by using kernel ridge regression.

We replicated the experiment one hundred times. In each replicate, we used stratified sampling to randomly pick $10\%$ of the observations as the testing set and the remaining as the training set. The Gaussian kernel function was used for the kernel ridge regression. We set the number of nodes $k=\{10,30,50,70,90,110\}$. We followed the procedure in [52] to select the bandwidth for the kernel and the regularization parameters. The performance of different estimators was evaluated by the MSE, respecting the prediction on the testing set. The experiment was conducted in R using a Windows computer with 16GB of memory and a single-threaded 3.5Ghz CPU.

The left panel of Fig. 6 shows the MSE versus the different number of nodes $k$, and the right panel of Fig. 6 shows the CPU time (in seconds) versus different $k$. To the best of our knowledge, our work is the first approach that addresses the response-skewed issue under the dacKRR setting, while other approaches only considered discrete-response settings or cannot be easily extended to the divide-and-conquer scenario. Thus we only compare our method with the classical dacKRR approach. In Fig. 6, the blue lines and the red lines represent the classical dacKRR approach and the proposed approach, respectively. Vertical bars represent the standard errors, which are obtained from one hundred replicates. The black lines represent the results of the full sample KRR estimate, and the gray dashed lines represent its standard errors. We can see that the proposed estimate consistently outperforms the classical dacKRR estimate. In particular, we observe the MSE of the proposed estimate is comparable with the MSE of the full sample KRR estimate. The MSE for the classical dacKRR estimate, however, almost doubles the MSE for the full sample KRR estimate when the number of nodes $k$ is greater than one hundred. The bottom panel of Fig. 6 shows the CPU time of the proposed approach is comparable with the CPU time of the classical divide-and-conquer approach. All these observations are consistent with the findings in the previous section. Such observations indicate the proposed approach outperforms the classical dacKRR approach for response-skewed data without requiring too much extra computing time.