As of the last couple decades, the overwhelming data accumulation in various domains such as medical research, road safety, bioinformatics, physics, genetics and social media, led scientists to rethink standard statistical learning methods. Beyond the computational barrier (in time and memory) induced by the data deluge and the economic cost of training models on large datasets, the overconsumption of energy resources also became a concerning topic for the scientific community. With this in mind, it is of interest to reduce the computational cost of learning methods, particularly regarding their training phase, which often involves minimizing an objective function and is more demanding than inference.

A first thought to limit the computational burden of training procedures, is to leverage more efficient optimization algorithms, such as active set strategies (24), stochastic gradient methods (3; 7) and random coordinate descents (39; 14), to cite only a few. However, these approaches still work on the whole training set, which may include either redundant or quite uninformative observations. Data set reduction methods can then be used in order to reduce the effective size of the dataset, while producing only informative datapoints. Sketching, for instance, is a way to reduce the dataset by generating new individuals, which are linear combinations of the original ones (30; 32; 47). Another approach is the coreset technique, which randomly pick points from de dataset in order to uniformly approximate the objective function (22; 29; 25; 35; 5). Nonuniform subsampling of points can also be implemented thanks to their leverage scores (28; 9; 49; 27), which are based on the covariates but not on the target values. Stepping back, subsampling can also be analyzed for general inclusion probabilities, from the excess risk point of view (8; 18) and the estimator asymptotic distribution side (44; 46; 19). In the latter case, it is even possible to exhibit an optimal subsampling scheme, which minimizes the trace of the limit covariance matrix. Several other approaches exist, such as sampling with influence functions (41) and leveraging minimum enclosing ball from computational geometry (42).

At this step, we remark that many methods, such as stochastic optimization and sketching, are task-agnostic, which leaves an area of improvement. Moreover, if coreset techniques propose data-dependent sampling procedures, these ones are aimed at uniform approximation the objective function and consequently often leverage intractable probabilities (29; 5). On the contrary, optimal subsampling proposes an easy-to-implement way to extract informative individuals from the original dataset, such that the trace of the asymptotic covariance matrix is minimal (46).

In this work, we focus on nonparametric supervised learning, and in particular on kernel methods, where the hypotheses lie in a reproducing kernel Hilbert space (RKHS) (2). In this case, few of the methods described above can be used (they have been mainly designed for parametric inference), the most notable being sketching (50; 12; 11) and specific approaches tailored for support vector machines, such as minimum enclosing balls (42) and stochastic subgradient algorithm (38). However, kernel learning can rely on kernel approximation techniques to lighten the numerical training cost, which consists in replacing either the Gram matrix by a low-rank approximation (for instance, thanks to the Nyström method (48; 10; 33)), or the feature map by a finite-dimensional mapping (thanks to random Fourier features (RFF) (31; 4; 34; 26)). Since the computational cost of a kernel method is mainly function of the size of the optimization problem to solve (which is basically represented by the spatial complexity of the design or Gram matrix of the effective dataset), the appeal of the Nyström method, RFF and sketching is to cut down the computational burden from $O(n^{2})$ (for an initial dataset of size $n$) to $O(nm)$, where $m$ an approximation parameter (the dimension of the subspace for the Nyström method, the number of random projections for RFF, the number of linear combinations for the sketching approach). If this is a big gain, we remark that the computational cost still depends on the size $n$ of the original dataset, which may be very large in practice. On the contrary, subsampling $s$ datapoints would lead to a $O(s^{2})$ cost, which is largely less.

We thus aim at designing a data-dependent subsampling procedure for nonparametric supervised learning in RKHSs. From the methodological side, we follow the line of research initiated by 44 (then followed by numerous papers including (45; 46; 19; 40)), which determines the subsampling scheme minimizing the trace of the asymptotic covariance matrix of a M-estimator of an empirical risk reweighted à la Horvitz-Thompson. From the theoretical side, we study asymptotics in RKHSs with tools coming from functional random probability spaces (6; 21) with pieces of elements from martingale theory (17). These two aspects distinguish our work from that by 46 since we place ourself in a particular nonparametric setting and employ different techniques of proof. The latter are also different from the empirical process techniques used by 16; 36 for analyzing usual empirical risk minimizers (ERMs).

In Section 2, we first revisit asymptotic normality of the ERM in RKHSs, based on the asymptotically finite-dimensional property (21; 43). This result corroborates the ones presented in (16; 36) and comes with a complementary technique of proof. Then, in Section 3, the subsampled risk minimizer (SRM) is presented as well as its convergence and asymptotic distribution. Section 4 exhibits the optimal subsampling scheme and Section 5 presents a way to implement it in practice. Le latter strategy is numerically evaluated on synthetic and real-world datasets in Section 6, and favorably compared to state-of-the-art approximation techniques for kernel methods.¹¹ 1 The Python source code for reproducing the numerical experiments is provided at https://gitlab.com/Yaloude/rkhs-subsampling. All proofs are gathered in Appendix A.

Let $Z\in\mathcal{Z}\subset\mathbb{R}^{d+1}$ be the random variable of interest with unknown distribution. We focus here on supervised learning and assume that $Z$ has a canonical decomposition: $Z=(X,Y)\in\mathcal{X}\times\mathcal{Y}\subset\mathbb{R}^{d}\times\mathbb{R}$ where $X$ is the vector of covariates and $Y$ is the response. Denote $\ell:\mathbb{R}\times\mathcal{Y}\rightarrow\mathbb{R}$ the loss function associated to the learning task of interest, and assume that $\ell$ is non-negative, convex and twice-differentiable in its first argument. The setting defined this way suits both binary classification: $\mathcal{Y}=\{\pm 1\}$ and for instance $\ell(\hat{y},y)=\log(1+\exp(-y\hat{y}))$ for logistic regression (where $\hat{y}\in\mathbb{R}$ is a prediction compared to $y\in\mathcal{Y}$); as well as regression: $\mathcal{Y}=\mathbb{R}$ and for example $\ell(\hat{y},y)=\frac{1}{p}|\hat{y}-y|^{p}$ for $L^{p}$ regression.

The working set of predictors is $\mathbb{H}$, the RKHS induced by some kernel $k:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$ (i.e. a symmetric and positive definite function). In other words, $\mathbb{H}$ is the closure of the span of the map by $k$ of elements of $\mathcal{X}$, i.e.

where $K_{x}=k(x,\cdot)$ for all $x\in\mathcal{X}$.

Given the previous elements, the population risk function of the learning task considered is:

which could be minimized so as to obtain an optimal predictor. However two problems arise: first, the distribution of $Z=(X,Y)$ is assumed unknown (such that the expectation in (1) cannot be computed); second, the current assumptions do not ensure existence of a minimizer of $L$ in $\mathbb{H}$. To solve those problems, we turn towards the empirical version of (1):

based on a dataset ${\mathcal{D}_{n}}=\{Z_{1},\dots,Z_{n}\}=\{(X_{1},Y_{1}),\dots,(X_{n},Y_{n})\}$ made of independent and identically distributed (i.i.d.) random variables sharing the same distribution as $Z=(X,Y)$, and add a ridge penalization tuned by some parameter ${\lambda_{0}}>0$ for the true risk and by ${\lambda_{n}}>0$ (a possibly random parameter) for the empirical one. The objective functions, namely the true and empirical regularized risks now read:

and

Let us remark that minimizing ${\mathcal{R}}$ (respectively ${\mathcal{R}}_{n}$) is equivalent to minimizing $L$ (respectively $L_{n}$) inside a ball of $\mathbb{H}$, the radius of which depending both on the penalization parameter ${\lambda_{0}}$ (respectively ${\lambda_{n}}$) and the data (see for example (36)). In addition, the framework above offers the chance to take ${\lambda_{n}}\neq{\lambda_{0}}$, which makes sens in practice since penalizing is also a way to prevent overfitting, on the condition that the penalization parameter ${\lambda_{n}}$ is tuned with respect to the sample size $n$ and the data (roughly speaking, the less the data, the stronger the penalization).

Now, by convexity of the loss $\ell$ (with respect to its first argument) and ${\lambda_{0}}>0$, the true regularized risk ${\mathcal{R}}$ is lower semi-continuous and strongly convex, and admits a unique minimizer in $\mathbb{H}$, denoted ${h_{0}}$ (6, Lemma 2.15). The same holds true for the empirical regularized risk ${\mathcal{R}}_{n}$ (since ${\lambda_{n}}>0$), the minimizer of which is denoted ${\hat{h}_{n}}\in\mathbb{H}$.

We now introduce, in this section, the two considered subsampling strategies (with and without replacement), along with the corresponding estimators, nicknamed SRM in both cases. We then prove asymptotic normality of both estimators with general subsampling schemes, and discuss in the next section the choice of the latter.

Let $\Delta_{n-1}=\{\bm{x}\in(0,1)^{n},\sum_{i=1}^{n}x_{i}=1\}$ be the set of probability vectors of size $n$, and consider ${\bm{\pi}}_{n}\in\Delta_{n-1}$ a subsampling probability (SSP) vector, which may ultimately depend on the dataset ${\mathcal{D}_{n}}$. A subsampling scheme of size $s_{n}\leq n$ consists in randomly picking points $Z_{i}$ in the dataset ${\mathcal{D}_{n}}$ with inclusion probabilities $\pi_{n,i}$. More formally, we consider, in the context of this paper, the two following subsampling strategies: one with replacement and another without.

Since the role of subsampling is to reduce the dataset size by randomly picking points, it is of interest to choose these few points in order to maximize the SRM estimator precision. Once the limit ratio $c=\lim_{n\to\infty}s_{n}/n$ has been chosen (see Theorems 3.2 and 3.3), it is about minimizing the asymptotic covariance operator. The most common way to do it is to reduce the covariance information to a scalar, and to minimize the latter. Several manners of reducing the covariance information exist and are well documented in the literature of optimal design (see for instance (23)). In this study, we follow the work by 46 and focus on L-optimality, which consists in computing the trace of the covariance, preconditioned by a well-chosen operator (here $(\nabla^{2}R(h_{0}))^{-1}$, which exists according to (16, Lemma A.5)). This choice is consistent with the dependence of the covariance operator to $\psi$ and is computationally advantageous.

In our infinite dimensional framework based on RKHSs, the trace of a self-adjoint bounded linear map does not differ much from the finite dimensional case. For such an operator $\Lambda:\mathbb{H}\to\mathbb{H}$, it can be defined as $\Tr(\Lambda)\overset{\mathrm{def}}{=}\sum_{j=1}^{\infty}\langle\Lambda e_{j},e_{j}\rangle_{\mathbb{H}}$, and, in the particular case where $\Lambda$ is the outer product of two vectors, i.e. $\Lambda=f\otimes h$ (with $f,h\in\mathbb{H}$), then $\Tr(\Lambda)=\langle f,h\rangle_{\mathbb{H}}$ (see (21) for more details).

Going back to Theorems 3.3 and 3.2, let $\Sigma^{*}(\psi)$ be either $\Sigma^{*}_{\mathrm{M}}(\psi)$ or $\Sigma^{*}_{\mathrm{P},c}(\psi)$, depending on the subsampling strategy. In both cases, L-optimality boils down to minimizing

with respect to $\psi$, where $C$ is a constant term. Proposition 4.1 gives a closed-form expression of such an L-optimal absolute interest map $\psi$.

If it is clear that, given our technical Assumptions 1, 2 and 3, the L-optimal absolute interest function $\psi^{\star}$ verifies Assumption 6, it is not true neither for Assumption 5, nor for the basic requirement of Poisson subsampling (i.e. that all SSPs are bounded by $1/s_{n}$). This is discussed in the two forthcoming remarks.

Let us now investigate some examples of supervised learning settings, which illustrate the shape of the L-optimal absolute interest function $\psi^{\star}$, revealed in Proposition 4.1. For this purpose and for simplicity, we consider a radial basis function kernel $k$, implying that for all $\bm{x}\in\mathcal{X}$, $k(\bm{x},\bm{x})$ is constant.

We address in this section the practical side of subsampling with the L-optimal absolute interest function $\psi^{\star}$, exhibited in Proposition 4.1. First, as it can be noticed, $\psi^{\star}$ makes use of the unknown regularized risk minimizer ${h_{0}}$. Nevertheless, the latter can be estimated by plug-in as soon as we are provided with prior information on ${h_{0}}$, which is assumed to be true in the form of a naive estimator ${\hat{h}_{k_{n}}}$ of ${h_{0}}$ (called a pilot estimator, and computed with $k_{n}$ points). Based on it, L-optimal SSPs can be estimated as $\tilde{\pi}_{n,i}\propto|\ell^{\prime}({\hat{h}_{k_{n}}}(X_{i}),Y_{i})|\sqrt{k(X_{i},X_{i})}$ for all $i\in\{1,\dots,n\}$. Second, as explained in the previous section, there is no guarantee that the map $\psi^{\star}$ leads to SSPs which verify Assumption 5. For this reason, it is proposed to smooth estimated SSPs as:

for all $i\in\{1,\dots,n\}$, where $\alpha\in(0,1]$ is a hyperparameter. As a last point, since SSPs should be less than $1/s_{n}$ for Poisson subsampling, one may apply a post hoc truncation, such as $\tilde{\bm{\pi}}_{n}^{\alpha}\wedge s_{n}^{-1}$. Everything together provides data driven SSPs satisfying all requirements and usable to build an SRM ${h_{s_{n}}^{*\alpha}}$, called an $\alpha$-piloted subsampling estimator.

Theorems 5.1 and 5.2 below state that the latter estimator is consistent to ${h_{0}}$ and asymptotically normal. In particular, it exhibits a covariance operator (which is well-defined and finite by the very definition of $\psi^{\star}$ and Assumptions 1, 2 and 3), which, compared to that of Theorems 3.2 and 3.3, is slightly impacted by smoothing the SSPs (tuned by the hyperparameter $\alpha$) but not by the use of a plug-in estimation of $\psi^{\star}$.