Year: 2020

Abstract:
We study the average CVloo stability of kernel ridge-less regression and derive corresponding risk bounds. We show that the interpolating solution with minimum norm minimizes a bound on CVloo stability, which in turn is controlled by the condition number of the empirical kernel matrix. The latter can be characterized in the asymptotic regime where both the dimension and cardinality of the data go to infinity. Under the assumption of random kernel matrices, the corresponding test error should be expected to follow a double descent curve.

heading: 1 INTRODUCTION
text: Statistical learning theory studies the learning properties of machine learning algorithms, and more fundamentally, the conditions under which learning from finite data is possible. In this context, classical learning theory focuses on the size of the hypothesis space in terms of different complexity measures, such as combinatorial dimensions, covering numbers and Rademacher/Gaussian complexities (Shalev-Shwartz & Ben-David, 2014; Boucheron et al., 2005). Another more recent approach is based on defining suitable notions of stability with respect to perturbation of the data (Bousquet & Elisseeff, 2001; Kutin & Niyogi, 2002). In this view, the continuity of the process that maps data to estimators is crucial, rather than the complexity of the hypothesis space. Different notions of stability can be considered, depending on the data perturbation and metric considered (Kutin & Niyogi, 2002). Interestingly, the stability and complexity approaches to characterizing the learnability of problems are not at odds with each other, and can be shown to be equivalent as shown in Poggio et al. (2004) and Shalev-Shwartz et al. (2010).
In modern machine learning overparameterized models, with a larger number of parameters than the size of the training data, have become common. The ability of these models to generalize is well explained by classical statistical learning theory as long as some form of regularization is used in the training process (Bühlmann & Van De Geer, 2011; Steinwart & Christmann, 2008). However, it was recently shown - first for deep networks (Zhang et al., 2017), and more recently for kernel methods (Belkin et al., 2019) - that learning is possible in the absence of regularization, i.e., when perfectly fitting/interpolating the data. Much recent work in statistical learning theory has tried to find theoretical ground for this empirical finding. Since learning using models that interpolate is not exclusive to deep neural networks, we study generalization in the presence of interpolation in the case of kernel methods. We study both linear and kernel least squares problems in this paper.
Our Contributions:
• We characterize the generalization properties of interpolating solutions for linear and kernel least squares problems using a stability approach. While the (uniform) stability properties of regularized kernel methods are well known (Bousquet & Elisseeff, 2001), we study interpolating solutions of the unregularized ("ridgeless") regression problems.
• We obtain an upper bound on the stability of interpolating solutions, and show that this upper bound is minimized by the minimum norm interpolating solution. This also means that among all interpolating solutions, the minimum norm solution has the best test error. In
particular, the same conclusion is also true for gradient descent, since it converges to the minimum norm solution in the setting we consider, see e.g. Rosasco & Villa (2015). • Our stability bounds show that the average stability of the minimum norm solution is
controlled by the condition number of the empirical kernel matrix. It is well known that the numerical stability of the least squares solution is governed by the condition number of the associated kernel matrix (see the discussion of why overparametrization is “good” in Poggio et al. (2019)). Our results show that the condition number also controls stability (and hence, test error) in a statistical sense.
Organization: In section 2, we introduce basic ideas in statistical learning and empirical risk minimization, as well as the notation used in the rest of the paper. In section 3, we briefly recall some definitions of stability. In section 4, we study the stability of interpolating solutions to kernel least squares and show that the minimum norm solutions minimize an upper bound on the stability. In section 5 we discuss our results in the context of recent work on high dimensional regression. We conclude in section 6.

heading: 2 STATISTICAL LEARNING AND EMPIRICAL RISK MINIMIZATION
text: We begin by recalling the basic ideas in statistical learning theory. In this setting, X is the space of features, Y is the space of targets or labels, and there is an unknown probability distribution µ on the product space Z = X × Y . In the following, we consider X = Rd and Y = R. The distribution µ is fixed but unknown, and we are given a training set S consisting of n samples (thus |S| = n) drawn i.i.d. from the probability distribution on Zn, S = (zi)ni=1 = (xi, yi) n i=1. Intuitively, the goal of supervised learning is to use the training set S to “learn” a function fS that evaluated at a new value xnew should predict the associated value of ynew, i.e. ynew ≈ fS(xnew). The loss is a function V : F × Z → [0,∞), where F is the space of measurable functions from X to Y , that measures how well a function performs on a data point. We define a hypothesis space H ⊆ F where algorithms search for solutions. With the above notation, the expected risk of f is defined as I[f ] = EzV (f, z) which is the expected loss on a new sample drawn according to the data distribution µ. In this setting, statistical learning can be seen as the problem of finding an approximate minimizer of the expected risk given a training set S. A classical approach to derive an approximate solution is empirical risk minimization (ERM) where we minimize the empirical risk IS [f ] = 1 n ∑n i=1 V (f, zi).
A natural error measure for our ERM solution fS is the expected excess risk ES [I[fS ]−minf∈H I[f ]]. Another common error measure is the expected generalization error/gap given by ES [I[fS ]− IS [fS ]]. These two error measures are closely related since, the expected excess risk is easily bounded by the expected generalization error (see Lemma 5).

heading: 2.1 KERNEL LEAST SQUARES AND MINIMUM NORM SOLUTION
text: The focus in this paper is on the kernel least squares problem. We assume the loss function V is the square loss, that is, V (f, z) = (y − f(x))2. The hypothesis space is assumed to be a reproducing kernel Hilbert space, defined by a positive definite kernel K : X ×X → R or an associated feature map Φ : X → H, such that K(x,x′) = 〈Φ(x),Φ(x′)〉H for all x,x′ ∈ X , where 〈·, ·〉H is the inner product in H. In this setting, functions are linearly parameterized, that is there exists w ∈ H such that f(x) = 〈w,Φ(x)〉H for all x ∈ X . The ERM problem typically has multiple solutions, one of which is the minimum norm solution:
f†S = arg min f∈M ‖f‖H , M = arg min f∈H
1
n n∑ i=1 (f(xi)− yi)2. (1)
Here ‖·‖H is the norm onH induced by the inner product. The minimum norm solution can be shown to be unique and satisfy a representer theorem, that is for all x ∈ X:
f†S(x) = n∑ i=1 K(x,xi)cS [i], cS = K †y (2)
where cS = (cS [1], . . . , cS [n]),y = (y1 . . . yn) ∈ Rn, K is the n by n matrix with entries Kij = K(xi,xj), i, j = 1, . . . , n, and K† is the Moore-Penrose pseudoinverse of K. If we assume n ≤ d and that we have n linearly independent data features, that is the rank of X is n, then it is possible to show that for many kernels one can replace K† by K−1 (see Remark 2). Note that invertibility is necessary and sufficient for interpolation. That is, if K is invertible, f†S(xi) = yi for all i = 1, . . . , n, in which case the training error in (1) is zero.
Remark 1 (Pseudoinverse for underdetermined linear systems) A simple yet relevant example are linear functions f(x) = w>x, that correspond toH = Rd and Φ the identity map. If the rank of X ∈ Rd×n is n, then any interpolating solution wS satisfies w>S xi = yi for all i = 1, . . . , n, and the minimum norm solution, also called Moore-Penrose solution, is given by (w†S)
> = y>X† where the pseudoinverse X† takes the form X† = X>(XX>)−1.
Remark 2 (Invertibility of translation invariant kernels) Translation invariant kernels are a family of kernel functions given by K(x1,x2) = k(x1 − x2) where k is an even function on Rd. Translation invariant kernels are Mercer kernels (positive semidefinite) if the Fourier transform of k(·) is non-negative. For Radial Basis Function kernels (K(x1,x2) = k(||x1 − x2||)) we have the additional property due to Theorem 2.3 of Micchelli (1986) that for distinct points x1,x2, . . . ,xn ∈ Rd the kernel matrix K is non-singular and thus invertible.
The above discussion is directly related to regularization approaches.
Remark 3 (Stability and Tikhonov regularization) Tikhonov regularization is used to prevent potential unstable behaviors. In the above setting, it corresponds to replacing Problem (1) by minf∈H 1 n ∑n i=1(f(xi) − yi)2 + λ ‖f‖ 2 H where the corresponding unique solution is given by
fλS (x) = ∑n i=1K(x,xi)c[i], c = (K + λIn)
−1y. In contrast to ERM solutions, the above approach prevents interpolation. The properties of the corresponding estimator are well known. In this paper, we complement these results focusing on the case λ→ 0.
Finally, we end by recalling the connection between minimum norm and the gradient descent.
Remark 4 (Minimum norm and gradient descent) In our setting, it is well known that both batch and stochastic gradient iterations converge exactly to the minimum norm solution when multiple solutions exist, see e.g. Rosasco & Villa (2015). Thus, a study of the properties of the minimum norm solution explains the properties of the solution to which gradient descent converges. In particular, when ERM has multiple interpolating solutions, gradient descent converges to a solution that minimizes a bound on stability, as we show in this paper.

heading: 3 ERROR BOUNDS VIA STABILITY
text: In this section, we recall basic results relating the learning and stability properties of Empirical Risk Minimization (ERM). Throughout the paper, we assume that ERM achieves a minimum, albeit the extension to almost minimizer is possible (Mukherjee et al., 2006) and important for exponential-type loss functions (Poggio, 2020). We do not assume the expected risk to achieve a minimum. Since we will be considering leave-one-out stability in this section, we look at solutions to ERM over the complete training set S = {z1, z2, . . . , zn} and the leave one out training set Si = {z1, z2, . . . , zi−1, zi+1, . . . , zn} The excess risk of ERM can be easily related to its stability properties. Here, we follow the definition laid out in Mukherjee et al. (2006) and say that an algorithm is Cross-Validation leave-one-out (CVloo) stable in expectation, if there exists βCV > 0 such that for all i = 1, . . . , n,
ES [V (fSi , zi)− V (fS , zi)] ≤ βCV . (3) This definition is justified by the following result that bounds the excess risk of a learning algorithm by its average CVloo stability (Shalev-Shwartz et al., 2010; Mukherjee et al., 2006).
Lemma 5 (Excess Risk & CVloo Stability) For all i = 1, . . . , n, ES [I[fSi ]− inf
f∈H I[f ]] ≤ ES [V (fSi , zi)− V (fS , zi)]. (4)
Remark 6 (Connection to uniform stability and other notions of stability) Uniform stability, introduced by Bousquet & Elisseeff (2001), corresponds in our notation to the assumption that there exists βu > 0 such that for all i = 1, . . . , n, supz∈Z |V (fSi , z)− V (fS , z)| ≤ βu. Clearly this is a strong notion implying most other definitions of stability. We note that there are number of different notions of stability. We refer the interested reader to Kutin & Niyogi (2002) , Mukherjee et al. (2006).
We recall the proof of Lemma 5 in Appendix A.2 due to lack of space. In Appendix A, we also discuss other definitions of stability and their connections to concepts in statistical learning theory like generalization and learnability.
4 CVloo STABILITY OF KERNEL LEAST SQUARES
In this section we analyze the expected CVloo stability of interpolating solutions to the kernel least squares problem, and obtain an upper bound on their stability. We show that this upper bound on the expected CVloo stability is smallest for the minimum norm interpolating solution (1) when compared to other interpolating solutions to the kernel least squares problem.
We have a dataset S = {(xi, yi)}ni=1 and we want to find a mapping f ∈ H, that minimizes the empirical least squares risk. Here H is a reproducing kernel hilbert space (RKHS) defined by a positive definite kernel K : X × X → R. All interpolating solutions are of the form f̂S(·) =∑n j=1 ĉS [j]K(xj , ·), where ĉS = K†y + (I −K†K)v. Similarly, all interpolating solutions on
the leave one out dataset Si can be written as f̂Si(·) = ∑n j=1,j 6=i ĉSi [j]K(xj , ·), where ĉSi = K†Siyi + (I−K † Si KSi)vi. Here K,KSi are the empirical kernel matrices on the original and leave one out datasets respectively. We note that when v = 0 and vi = 0, we obtain the minimum norm interpolating solutions on the datasets S and Si.
Theorem 7 (Main Theorem) Consider the kernel least squares problem with a bounded kernel and bounded outputs y, that is there exist κ,M > 0 such that
K(x,x′) ≤ κ2, |y| ≤M, (5)
almost surely. Then for any interpolating solutions f̂Si , f̂S ,
ES [V (f̂Si , zi)− V (f̂S , zi)] ≤ βCV (K†,y,v,vi) (6) This bound βCV is minimized when v = vi = 0, which corresponds to the minimum norm interpolating solutions f†S , f † Si . For the minimum norm solutions we have
βCV = C1β1 + C2β2, where β1 = ES [ ||K 12 ||op||K†||op × cond(K)× ||y|| ] and, β2 =
ES [ ||K 12 ||2op||K†||2op × (cond(K))2 × ||y||2 ] , andC1, C2 are absolute constants that do not depend
on either d or n.
In the above theorem ||K||op refers to the operator norm of the kernel matrix K, ||y|| refers to the standard `2 norm for y ∈ Rn, and cond(K) is the condition number of the matrix K. We can combine the above result with Lemma 5 to obtain the following bound on excess risk for minimum norm interpolating solutions to the kernel least squares problem:
Corollary 8 The excess risk of the minimum norm interpolating kernel least squares solution can be bounded as: ES [ I[f†Si ]− inff∈H I[f ] ] ≤ C1β1 + C2β2
where β1, β2 are as defined previously.
Remark 9 (Underdetermined Linear Regression) In the case of underdetermined linear regression, ie, linear regression where the dimensionality is larger than the number of samples in the training set, we can prove a version of Theorem 7 with β1 = ES [∥∥X†∥∥ op ‖y‖ ] and
β2 = ES [∥∥X†∥∥2 op ‖y‖2 ] . Due to space constraints, we present the proof of the results in the linear regression case in Appendix B.

heading: 4.1 KEY LEMMAS
text: In order to prove Theorem 7 we make use of the following lemmas to bound the CVloo stability using the norms and the difference of the solutions.
Lemma 10 Under assumption (5), for all i = 1. . . . , n, it holds that ES [V (f̂Si , zi)− V (f̂S , zi)] ≤ ES [( 2M + κ (∥∥∥f̂S∥∥∥ H + ∥∥∥f̂Si∥∥∥H))× κ∥∥∥f̂S − f̂Si∥∥∥H]
Proof We begin, recalling that the square loss is locally Lipschitz, that is for all y, a, a′ ∈ R, with |(y − a)2 − (y − a′)2| ≤ (2|y|+ |a|+ |a′|))|a− a′|.
If we apply this result to f, f ′ in a RKHSH, |(y − f(x))2 − (y − f ′(x))2| ≤ κ(2M + κ (‖f‖H + ‖f ′‖H)) ‖f − f ′‖H .
using the basic properties of a RKHS that for all f ∈ H |f(x)| ≤ ‖f‖∞ = supx|f(x)| = supx|〈f,Kx〉H| ≤ κ ‖f‖H (7)
In particular, we can plug f̂Si and f̂S into the above inequality, and the almost positivity of ERM (Mukherjee et al., 2006) will allow us to drop the absolute value on the left hand side. Finally the desired result follows by taking the expectation over S.
Now that we have bounded the CVloo stability using the norms and the difference of the solutions, we can find a bound on the difference between the solutions to the kernel least squares problem. This is our main stability estimate.
Lemma 11 Let f̂S , f̂Si be any interpolating kernel least squares solutions on the full and leave one out datasets (as defined at the top of this section), then ∥∥∥f̂S − f̂Si∥∥∥H ≤ BCV (K†,y,v,vi), and BCV is minimized when v = vi = 0, which corresponds to the minimum norm interpolating solutions f†S , f † Si . Also for some absolute constant C,∥∥∥f†S − f†Si∥∥∥H ≤ C × ∥∥∥K 12 ∥∥∥op ∥∥K†∥∥op × cond(K)× ‖y‖ (8) Since the minimum norm interpolating solutions minimize both ∥∥∥f̂S∥∥∥ H + ∥∥∥f̂Si∥∥∥H and ∥∥∥f̂S − f̂Si∥∥∥H (from lemmas 10, 11), we can put them together to prove theorem 7. In the following section we provide the proof of Lemma 11.
Remark 12 (Zero training loss) In Lemma 10 we use the locally Lipschitz property of the squared loss function to bound the leave one out stability in terms of the difference between the norms of the solutions. Under interpolating conditions, if we set the term V (f̂S , zi) = 0, the leave one
out stability reduces to ES [ V (f̂Si , zi)− V (f̂S , zi) ] = ES [ V (f̂Si , zi) ] = ES [(f̂Si(xi)− yi)2] =
ES [(f̂Si(xi)− f̂S(xi))2] = ES [〈f̂Si(·)− f̂S(·),Kxi(·)〉2] ≤ ES [ ||f̂S − f̂Si ||2H × κ2 ] . We can plug
in the bound from Lemma 11 to obtain similar qualitative and quantitative (up to constant factors) results as in Theorem 7.
Simulation: In order to illustrate that the minimum norm interpolating solution is the best performing interpolating solution we ran a simple experiment on a linear regression problem. We synthetically generated data from a linear model y = w>X, where X ∈ Rd×n was i.i.d N (0, 1). The dimension of the data was d = 1000 and there were n = 200 samples in the training dataset. A held out dataset of 50 samples was used to compute the test mean squared error (MSE). Interpolating solutions were computed as ŵ> = y>X†+v>(I−XX†) and the norm of v was varied to obtain the plot. The results are shown in Figure 1, where we can see that the training loss is 0 for all interpolants, but test MSE increases as ||v|| increases, with (w†)> = y>X† having the best performance. The figure reports results averaged over 100 trials.

heading: 4.2 PROOF OF LEMMA 11
text: We can write any interpolating solution to the kernel regression problem as f̂S(x) =∑n i=1 ĉS [i]K(xi,x) where ĉS = K
†y + (I − K†K)v, and K ∈ Rn×n is the kernel matrix K on S and v is any vector in Rn. i.e. Kij = K(xi,xj), and y ∈ Rn is the vector y = [y1 . . . yn]>. Similarly, the coefficient vector for the corresponding interpolating solution to the problem over the leave one out dataset Si is ĉSi = (KSi)
†yi + (I− (KSi)†KSi)vi. Where yi = [y1, . . . , 0, . . . yn]> and KSi is the kernel matrix K with the i
th row and column set to zero, which is the kernel matrix for the leave one out training set.
We define a = [−K(x1,xi), . . . ,−K(xn,xi)]> ∈ Rn and b ∈ Rn as a one-hot column vector with all zeros apart from the ith component which is 1. Let a∗ = a +K(xi,xi)b. Then, we have:
K∗ = K + ba > ∗
KSi = K∗ + ab > (9)
That is, we can write KSi as a rank-2 update to K. This can be verified by simple algebra, and using the fact that K is a symmetric kernel. Now we are interested in bounding ||f̂S − f̂Si ||H. For a function h(·) = ∑m i=1 piK(xi, ·) ∈ H we have ||h||H = √ p>Kp = ||K 12p||. So we have:
||f̂S − f̂Si ||H = ||K 1 2 (ĉS − ĉSi)||
= ||K 12 (K†y + (I−K†K)v − (KSi)†yi − (I− (KSi)†KSi)vi)||
= ||K 12 (K†y − (KSi)†y + yi(KSi)†b + (I−K†K)(v − vi) + (K†K− (KSi)†KSi)vi)||
= ||K 12 [(K† − (KSi)†)y + (I−K†K)(v − vi)− (K†K− (KSi)†KSi)vi]|| (10)
Here we make use of the fact that (KSi) †b = 0. If K has full rank (as in Remark 2), we see that b lies in the column space of K and a∗ lies in the column space of K>. Furthermore, β∗ = 1 +a>∗K †b = 1 +a>K†b+K(xi,xi)b >K†b = Kii(K †)ii 6= 0. Using equation 2.2 of Baksalary
et al. (2003) we obtain:
K†∗ = K † − (Kii(K†)ii)−1K†ba>∗K†
= K† − (Kii(K†)ii)−1K†ba>K† − ((K†)ii)−1K†bb>K†
= K† + (Kii(K †)ii) −1K†bb> − ((K†)ii)−1K†bb>K† (11)
Here we make use of the fact that a>K† = −b. Also, using the corresponding formula from List 2 of Baksalary et al. (2003), we have K†∗K∗ = K†K.
Next, we see that since K∗ has the same rank as K, a lies in the column space of K∗, and b lies in the column space of K>∗ . Furthermore β = 1 + b
>K∗a = 0. This means we can use Theorem 6 in Meyer (1973) (equivalent to formula 2.1 in Baksalary et al. (2003)) to obtain the expression for (KSi) †, with k = K†∗a and h = b>K † ∗.
(KSi) † = K†∗ − kk†K†∗ −K†∗h†h + (k†K†∗h†)kh
=⇒ (KSi)† −K†∗ = (k†K†∗h†)kh− kk†K†∗ −K†∗h†h =⇒ ||(KSi)† −K†∗||op ≤ 3||K†∗||op
(12)
Above, we use the fact that the operator norm of a rank 1 matrix is given by ||uv>||op = ||u|| × ||v||. Also, using the corresponding formula from List 2 of Baksalary et al. (2003), we have:
(KSi) †KSi = K † ∗K∗ − kk†
=⇒ K†K− (KSi)†KSi = kk† (13)
Putting the two parts together we obtain the bound on ∥∥(KSi)† −K†∥∥op:
||K† − (KSi)†||op = ||K† −K†∗ + K†∗ − (KSi)†||op ≤ 3||K†∗||op + ||K† −K†∗||op ≤ 3||K†||op + 4(Kii(K†)ii)−1||K†||op + 4((K†)ii)−1||K†||2op ≤ ||K†||op(3 + 8||K†||op||K||op)
(14)
The last step follows from (Kii)−1 ≤ ||K†||op and ((K†)ii)−1 ≤ ||K||op. Plugging in these calculations into equation 10 we get:
||f̂S − f̂Si ||H = ||K 1 2 [(K† − (KSi)†)y + (I−K†K)(v − vi)− (K†K− (KSi)†KSi)vi]|| ≤ ||K 12 ||op ( ||(K† − (KSi)†)y||+ ||(I−K†K)(v − vi)||+ ||kk†vi|| ) ≤ ||K 12 ||op(B0 + ||I−K†K||op||v − vi||+ ||vi||)
(15)
We see that the right hand side is minimized when v = vi = 0. We have also computed B0 = C × ||K†||op × cond(K)× ||y||, which concludes the proof of Lemma 11.

heading: 5 REMARK AND RELATED WORK
text: In the previous section we obtained bounds on the CVloo stability of interpolating solutions to the kernel least squares problem. Our kernel least squares results can be compared with stability bounds for regularized ERM (see Remark 3). Regularized ERM has a strong stability guarantee in terms of a uniform stability bound which turns out to be inversely proportional to the regularization parameter λ and the sample size n (Bousquet & Elisseeff, 2001). However, this estimate becomes vacuous as λ→ 0. In this paper, we establish a bound on average stability, and show that this bound is minimized when the minimum norm ERM solution is chosen. We study average stability since one can expect
20 40 60 80 100 120 140
n
0
500
1000
1500
2000
2500
3000
3500
C on
di tio
n nu
m be
r
d=nd>n d<n
Figure 2: Typical double descent of the condition number (y axis) of a radial basis function kernel K(x, x′) = exp ( − ||x−x ′||2
2σ2
) built from a random data matrix distributed asN (0, 1): as in the linear
case, the condition number is worse when n = d, better if n > d (on the right of n = d) and also better if n < d (on the left of n = d). The parameter σ was chosen to be 5. From Poggio et al. (2019)
worst case scenarios where the minimum norm is arbitrarily large (when n ≈ d). One of our key findings is the relationship between minimizing the norm of the ERM solution and minimizing a bound on stability.
This leads to a second observation, namely, that we can consider the limit of our risk bounds as both the sample size (n) and the dimensionality of the data (d) go to infinity, but the ratio dn → γ > 1 as n, d→∞ . This is a classical setting in statistics which allows us to use results from random matrix theory (Marchenko & Pastur, 1967). In particular, for linear kernels the behavior of the smallest eigenvalue of the kernel matrix (which appears in our bounds) can be characterized in this asymptotic limit. In fact, under appropriate distributional assumptions, our bound for linear regression can be computed as (||X†|| × ||y||)2 ≈ √ n√
d− √ n → 1√γ−1 . Here the dimension of the data coincides with
the number of parameters in the model. Interestingly, analogous results hold for more general kernels (inner product and RBF kernels) (El Karoui, 2010) where the asymptotics are taken with respect to the number and dimensionality of the data. These results predict a double descent curve for the condition number as found in practice, see Figure 2. While it may seem that our bounds in Theorem 7 diverge if d is held constant and n→∞, this case is not covered by our theorem, since when n > d we no longer have interpolating solutions.
Recently, there has been a surge of interest in studying linear and kernel least squares models, since classical results focus on situations where constraints or penalties that prevent interpolation are added to the empirical risk. For example, high dimensional linear regression is considered in Mei & Montanari (2019); Hastie et al. (2019); Bartlett et al. (2019), and “ridgeless” kernel least squares is studied in Liang et al. (2019); Rakhlin & Zhai (2018) and Liang et al. (2020). While these papers study upper and lower bounds on the risk of interpolating solutions to the linear and kernel least squares problem, ours are the first to derived using stability arguments. While it might be possible to obtain tighter excess risk bounds through careful analysis of the minimum norm interpolant, our simple approach helps us establish a link between stability in statistical and in numerical sense.
Finally, we can compare our results with observations made in Poggio et al. (2019) on the condition number of random kernel matrices. The condition number of the empirical kernel matrix is known to control the numerical stability of the solution to a kernel least squares problem. Our results show that the statistical stability is also controlled by the condition number of the kernel matrix, providing a natural link between numerical and statistical stability.

heading: 6 CONCLUSIONS
text: In summary, minimizing a bound on cross validation stability minimizes the expected error in both the classical and the modern regime of ERM. In the classical regime (d < n), CVloo stability implies generalization and consistency for n→∞. In the modern regime (d > n), as described in this paper, CVloo stability can account for the double descent curve in kernel interpolants (Belkin et al., 2019) under appropriate distributional assumptions. The main contribution of this paper is characterizing stability of interpolating solutions, in particular deriving excess risk bounds via a stability argument. In the process, we show that among all the interpolating solutions, the one with minimum norm also minimizes a bound on stability. Since the excess risk bounds of the minimum norm interpolant depend on the pseudoinverse of the kernel matrix, we establish here an elegant link between numerical and statistical stability. This also holds for solutions computed by gradient descent, since gradient descent converges to minimum norm solutions in the case of “linear” kernel methods. Our approach is simple and combines basic stability results with matrix inequalities.

heading: A EXCESS RISK, GENERALIZATION, AND STABILITY
text: We use the same notation as introduced in Section 2 for the various quantities considered in this section. That is in the supervised learning setup V (f, z) is the loss incurred by hypothesis f on the sample z, and I[f ] = Ez[V (f, z)] is the expected error of hypothesis f . Since we are interested in different forms of stability, we will consider learning problems over the original training set S = {z1, z2, . . . , zn}, the leave one out training set Si = {z1, . . . , zi−1, zi+1, . . . , zn}, and the replace one training set (Si, z) = {z1, . . . , zi−1, zi+1, . . . , zn, z}
A.1 REPLACE ONE AND LEAVE ONE OUT ALGORITHMIC STABILITY
Similar to the definition of expected CVloo stability in equation (3) of the main paper, we say an algorithm is cross validation replace one stable (in expectation), denoted as CVro, if there exists βro > 0 such that
ES,z[V (fS , z)− V (f(Si,z), z)] ≤ βro.
We can strengthen the above stability definition by introducing the notion of replace one algorithmic stability (in expectation) Bousquet & Elisseeff (2001). There exists αro > such that for all i = 1, . . . , n,
ES,z[ ∥∥fS − f(Si,z)∥∥∞] ≤ αro.
We make two observations: First, if the loss is Lipschitz, that is if there exists CV > 0 such that for all f, f ′ ∈ H
‖V (f, z)− V (f ′, z)‖ ≤ CV ‖f − f ′‖ ,
then replace one algorithmic stability implies CVro stability with βro = CV αro. Moreover, the same result holds if the loss is locally Lipschitz and there exists R > 0, such that ‖fS‖∞ ≤ R almost surely. In this latter case the Lipschitz constant will depend on R. Later, we illustrate this situation for the square loss.
Second, we have for all i = 1, . . . , n, S and z, ES,z[ ∥∥fS − f(Si,z)∥∥∞] ≤ ES,z[‖fS − fSi‖∞] + ES,z[∥∥f(Si,z) − fSi∥∥∞].
This observation motivates the notion of leave one out algorithmic stability (in expectation) Bousquet & Elisseeff (2001)]
ES,z[‖fS − fSi‖∞] ≤ αloo.
Clearly, leave one out algorithmic stability implies replace one algorithmic stability with αro = 2αloo and it implies also CVro stability with βro = 2CV αloo.
A.2 EXCESS RISK AND CVloo, CVro STABILITY
We recall the statement of Lemma 5 in section 3 that bounds the excess risk using the CVloo stability of a solution.
Lemma 13 (Excess Risk & CVloo Stability) For all i = 1, . . . , n,
ES [I[fSi ]− inf f∈H I[f ]] ≤ ES [V (fSi , zi)− V (fS , zi)]. (16)
In this section, two properties of ERM are useful, namely symmetry, and a form of unbiasedeness.
Symmetry. A key property of ERM is that it is symmetric with respect to the data set S, meaning that it does not depend on the order of the data in S.
A second property relates the expected ERM with the minimum of expected risk.
ERM Bias. The following inequality holds.
E[[IS [fS ]]−min f∈H I[f ] ≤ 0. (17)
To see this, note that IS [fS ] ≤ IS [f ]
for all f ∈ H by definition of ERM, so that taking the expectation of both sides ES [IS [fS ]] ≤ ES [IS [f ]] = I[f ]
for all f ∈ H. This implies ES [IS [fS ]] ≤ min
f∈H I[f ]
and hence (17) holds.
Remark 14 Note that the same argument gives more generally that
E[ inf f∈H [IS [f ]]− inf f∈H I[f ] ≤ 0. (18)
Given the above premise, the proof of Lemma 5 is simple.
Proof [of Lemma 5] Adding and subtracting ES [IS [fS ]] from the expected excess risk we have that ES [I[fSi ]−min
f∈H I[f ]] = ES [I[fSi ]− IS [fS ] + IS [fS ]−min f∈H I[f ]], (19)
and since ES [IS [fS ]]−minf∈H I[f ]] is less or equal than zero, see (18), then ES [I[fSi ]−min
f∈H I[f ]] ≤ ES [I[fSi ]− IS [fS ]]. (20)
Moreover, for all i = 1, . . . , n
ES [I[fSi ]] = ES [EziV (fSi , zi)] = ES [V (fSi , zi)] and
ES [IS [fS ]] = 1
n n∑ i=1 ES [V (fS , zi)] = ES [V (fS , zi)].
Plugging these last two expressions in (20) and in (19) leads to (4).
We can prove a similar result relating excess risk with CVro stability.
Lemma 15 (Excess Risk & CVro Stability) Given the above definitions, the following inequality holds for all i = 1, . . . , n,
ES [I[fS ]− inf f∈H I[f ]] ≤ ES [I[fS ]− IS [fS ]] = ES,z[V (fS , z)− V (f(Si,z), z)]. (21)
Proof The first inequality is clear from adding and subtracting IS [fS ] from the expected risk I[fS ] we have that
ES [I[fS ]−min f∈H I[f ]] = ES [I[fS ]− IS [fS ] + IS [fS ]−min f∈H I[f ]],
and recalling (18). The main step in the proof is showing that for all i = 1, . . . , n,
E[IS [fS ]] = E[V (f(Si,z), z)] (22)
to be compared with the trivial equality, E[IS [fS ] = E[V (fS , zi)]. To prove Equation (22), we have for all i = 1, . . . , n,
ES [IS [fS ]] = ES,z[ 1
n n∑ i=1 V (fS , zi)] = 1 n n∑ i=1 ES,z[V (f(Si,z), z)] = ES,z[V (f(Si,z), z)]
where we used the fact that by the symmetry of the algorithm ES,z[V (f(Si,z), z)] is the same for all i = 1, . . . , n. The proof is concluded noting that ES [I[fS ]] = ES,z[V (fS , z)].
A.3 DISCUSSION ON STABILITY AND GENERALIZATION
Below we discuss some more aspects of stability and its connection to other quantities in statistical learning theory.
Remark 16 (CVloo stability in expectation and in probability) In Mukherjee et al. (2006), CVloo stability is defined in probability, that is there exists βPCV > 0, 0 < δ P CV ≤ 1 such that
PS{|V (fSi , zi)− V (fS , zi)| ≥ βPCV } ≤ δPCV .
Note that the absolute value is not needed for ERM since almost positivity holds Mukherjee et al. (2006), that is V (fSi , zi)− V (fS , zi) > 0. Then CVloo stability in probability and in expectation are clearly related and indeed equivalent for bounded loss functions. CVloo stability in expectation (3) is what we study in the following sections.
Remark 17 (Connection to uniform stability and other notions of stability) Uniform stability, introduced by Bousquet & Elisseeff (2001), corresponds in our notation to the assumption that there exists βu > 0 such that for all i = 1, . . . , n, supz∈Z |V (fSi , z) − V (fS , z)| ≤ βu. Clearly this is a strong notion implying most other definitions of stability. We note that there are number of different notions of stability. We refer the interested reader to Kutin & Niyogi (2002) , Mukherjee et al. (2006).
Remark 18 (CVloo Stability & Learnability) A natural question is to which extent suitable notions of stability are not only sufficient but also necessary for controlling the excess risk of ERM. Classically, the latter is characterized in terms of a uniform version of the law of large numbers, which itself can be characterized in terms of suitable complexity measures of the hypothesis class. Uniform stability is too strong to characterize consistency while CVloo stability turns out to provide a suitably weak definition as shown in Mukherjee et al. (2006), see also Kutin & Niyogi (2002), Mukherjee et al. (2006). Indeed, a main result in Mukherjee et al. (2006) shows that CVloo stability is equivalent to consistency of ERM:
Theorem 19 Mukherjee et al. (2006) For ERM and bounded loss functions, CVloo stability in probability with βPCV converging to zero for n→∞ is equivalent to consistency and generalization of ERM.
Remark 20 (CVloo stability & in-sample/out-of-sample error) Let (S, z) = {z1, . . . , zn, z}, (z is a data point drawn according to the same distribution) and the corresponding ERM solution f(S,z), then (4) can be equivalently written as,
ES [I[fS ]− inf f∈F I[f ]] ≤ ES,z[V (fS , z)− V (f(S,z), z)].
Thus CVloo stability measures how much the loss changes when we test on a point that is present in the training set and absent from it. In this view, it can be seen as an average measure of the difference between in-sample and out-of-sample error.
Remark 21 (CVloo stability and generalization) A common error measure is the (expected) generalization gap ES [I[fS ]−IS [fS ]]. For non-ERM algorithms, CVloo stability by itself not sufficient to control this term, and further conditions are needed Mukherjee et al. (2006), since
ES [I[fS ]− IS [fS ]] = ES [I[fS ]− IS [fSi ]] + ES [IS [fSi ]− IS [fS ]].
The second term becomes for all i = 1, . . . , n,
ES [IS [fSi ]− IS [fS ]] = 1
n n∑ i=1 ES [V (fSi , zi)− V (fS , zi)] = ES [V (fSi , zi)− V (fS , zi)]
and hence is controlled by CV stability. The first term is called expected leave one out error in Mukherjee et al. (2006) and is controlled in ERM as n→∞, see Theorem 19 above.
B CVloo STABILITY OF LINEAR REGRESSION
We have a dataset S = {(xi, yi)}ni=1 and we want to find a mapping w ∈ R d, that minimizes the empirical least squares risk. All interpolating solutions are of the form ŵS = y>X†+v>(I−XX†). Similarly, all interpolating solutions on the leave one out dataset Si can be written as ŵSi = y>i (Xi) † + v>i (I−Xi(Xi)†). Here X,Xi ∈ R d×n are the data matrices for the original and leave one out datasets respectively. We note that when v = 0 and vi = 0, we obtain the minimum norm interpolating solutions on the datasets S and Si.
In this section we want to estimate the CVloo stability of the minimum norm solution to the ERM problem in the linear regression case. This is the case outlined in Remark 9 of the main paper. In order to prove Remark 9, we only need to combine Lemma 10 with the linear regression analogue of Lemma 11. We state and prove that result in this section. This result predicts a double descent curve for the norm of the pseudoinverse as found in practice, see Figure 3.
Lemma 22 Let ŵS , ŵSi be any interpolating least squares solutions on the full and leave one out datasets S, Si, then ‖ŵS − ŵSi‖ ≤ BCV (X†,y,v,vi), and BCV is minimized when v = vi = 0, which corresponds to the minimum norm interpolating solutions w†S ,w † Si . Also, ∥∥∥w†S −w†Si∥∥∥ ≤ 3∥∥X†∥∥op × ‖y‖ (23) As mentioned before in section 2.1 of the main paper, linear regression can be viewed as a case of the kernel regression problem whereH = Rd, and the feature map Φ is the identity map. The inner product and norms considered in this case are also the usual Euclidean inner product and 2-norm for vectors in Rd. The notation ‖·‖ denotes the Euclidean norm for vectors both in Rd and Rn. The usage of the norm should be clear from the context. Also, ‖A‖op is the left operator norm for a matrix A ∈ Rn×d, that is ‖A‖op = supy∈Rn,||y||=1 ||y>A||.
We have n samples in the training set for a linear regression problem, {(xi, yi)}ni=1. We collect all the samples into a single matrix/vector X = [x1x2 . . .xn] ∈ Rd×n, and y = [y1y2 . . . yn]> ∈ Rn. Then any interpolating ERM solution wS satisfies the linear equation
w>SX = y > (24)
Any interpolating solution can be written as:
(ŵS) > = y>X† + v>(I−XX†). (25)
If we consider the leave one out training set Si we can find the minimum norm ERM solution for Xi = [x1 . . .0 . . .xn] and yi = [y1 . . . 0 . . . yn]> as
(ŵSi) > = y>i (Xi) † + v>i (I−Xi(Xi)†). (26) We can write Xi as:
Xi = X + ab > (27)
where a ∈ Rd is a column vector representing the additive change to the ith column, i.e, a = −xi, and b ∈ Rn×1 is the i−th element of the canonical basis in Rn (all the coefficients are zero but the i−th which is one). Thus ab> is a d × n matrix composed of all zeros apart from the ith column which is equal to a.
We also have yi = y − yib. Now per Lemma 10 we are interested in bounding the quantity ||ŵSi − ŵS || = ||(ŵSi)> − (ŵS)>||. This simplifies to:
||ŵSi − ŵS || = ||y>i (Xi)† − y>X† + v>i − v> + v>XX† − v>i Xi(Xi)†|| = ||(y> − yib>)(Xi)† − y>X† + v>i − v> + v>XX† − v>i Xi(Xi)†|| = ||y>((Xi)† −X†) + yib>(Xi)† + v>i − v> + v>XX† − v>i Xi(Xi)†|| = ||y>((Xi)† −X†) + v>i − v> + v>XX† − v>i Xi(Xi)†|| = ||y>((Xi)† −X†) + (v>i − v>)(I−XX†)− v>i (XX† −Xi(Xi)†)|| (28)
In the above equation we make use of the fact that b>(Xi)† = 0. We use an old formula (Meyer, 1973; Baksalary et al., 2003) to compute (Xi)† from X†. We use the development of pseudo-inverses of perturbed matrices in Meyer (1973). We see that a = −xi is a vector in the column space of X and b is in the range space of XT (provided X has full column rank), with β = 1 + b>X†a = 1− b>X†xi = 0. This means we can use Theorem 6 in Meyer (1973) (equivalent to formula 2.1 in Baksalary et al. (2003)) to obtain the expression for (Xi)†
(Xi) † = X† − kk†X† −X†h†h + (k†X†h†)kh (29)
where k = X†a, and h = b>X†, and u† = u >
||u||2 for any non-zero vector u.
(Xi) † −X† = (k†X†h†)kh− kk†X† −X†h†h
= a>(X†)>X†(X†)>b× kh ||k||2||h||2 − kk†X† −X†h†h
=⇒ ||(Xi)† −X†||op ≤ |a>(X†)>X†(X†)>b| ||X†a||||b>X†|| + 2||X†||op
≤ ||X †||op||X†a||||b>X†|| ||X†a||||b>X†|| + 2||X†||op
= 3||X†||op
(30)
The above set of inequalities follows from the fact that the operator norm of a rank 1 matrix is given by ||uv>||op = ||u|| × ||v||
Also, from List 2 of Baksalary et al. (2003) we have that Xi(Xi)† = XX† − h†h. Plugging in these calculations into equation 28 we get:
||ŵSi − ŵS || = ||y>((Xi)† −X†) + (v>i − v>)(I−XX†)− v>i (XX† −Xi(Xi)†)|| ≤ B0 + ||I−XX†||op||v − vi||+ ||vi|| × ||h†h||op ≤ B0 + 2||v − vi||+ ||vi|| (31)
We see that the right hand side is minimized when v = vi = 0. We can also compute B0 = 3||X†||op||y||, which concludes the proof of Lemma 22.

Review summaries:
This paper investigates kernel ridge-less regression from a stability viewpoint by deriving its risk bounds. Using stability arguments to derive risk bounds have been widely adopting in machine learning. However, related studies on kernel ridge-less regression are still sparse. The present study fills this gap, which, in my opinion, is also one of the main contributions of the present study.