Implicit Regularization of Random Feature Models

Introduction

In this paper, we consider the Random Feature (RF) model which is an approximation of Kernel Methods (Rahimi & Recht, 2008) which has seen many recent theoretical developements.

The conventional wisdom suggests that to ensure good generalization performance, one should choose a model class that is complex enough to learn the signal from the training data, yet simple enough to avoid fitting spurious patterns therein (Bishop, 2006). This view has been questioned by recent developments in machine learning. First, Zhang et al. (2016) observed that modern neural network models can perfectly fit randomly labeled training data, while still generalizing well. Second, the test error as a function of parameters exhibits a so-called ‘double-descent’ curve for many models including neural networks, random forests, and random feature models (Advani & Saxe, 2017; Spigler et al., 2018; Belkin et al., 2018; Mei & Montanari, 2019; Belkin et al., 2019; Nakkiran et al., 2019).

The above models share the feature that for fixed input, the learned predictor $\hat{f}$ is random: for neural networks, this is due to the random initialization of the parameters and/or to the stochasticity of the training algorithm; for random forests, to the random branching; for random feature models, to the sampling of random features. The somehow surprising generalization behavior of these models has recently been the subject of increasing attention. In general, the risk (i.e. test error) is a random variable with two sources of randomness: the usual one due to the sampling of the training set, and the second one due to the randomness of the model itself.

We consider the Random Feature (RF) model (Rahimi & Recht, 2008) with features sampled from a Gaussian Process (GP) and study the RF predictor $\hat{f}$ minimizing the regularized least squares error, isolating the randomness of the model by considering fixed training data points. RF models have been the subject of intense research activity: they are (randomized) approximations of Kernel Methods aimed at easing the computational challenges of Kernel Methods while being asymptotically equivalent to them (Rahimi & Recht, 2008; Yang et al., 2012; Sriperumbudur & Szabó, 2015; Yu et al., 2016). Unlike the asymptotic behavior, which is well studied, RF models with a finite number of features are much less understood.

We consider a model of Random Features (RF) approximating a kernel method with kernel $K$. This model consists of $P$ Gaussian features, sampled i.i.d. from a (centered) Gaussian process with covariance kernel $K$. For a given training set of size $N$, we study the distribution of the RF predictor $\hat{f}^{(RF)}_{\lambda}$ with ridge parameter $\lambda>0$ ($L^{2}$ penalty on the parameters) and denote it by $\lambda$-RF. We show the following:

The distribution of $\hat{f}^{(RF)}_{\lambda}$ is that of a mixture of Gaussian processes.

2 Related works

Generalization of Random Features. The generalization behavior of Random Feature models has seen intense study in the Statistical Learning Theory framework. Rahimi & Recht (2009) find that $\mathcal{O}(N)$ features are sufficient to ensure the $\mathcal{O}(\frac{1}{\sqrt{N}})$ decay of the generalization error of Kernel Ridge Regression (KRR). Rudi & Rosasco (2017) improve on their result and show that $\mathcal{O}(\sqrt{N}\log N)$ features is actually enough to obtain the $\mathcal{O}(\frac{1}{\sqrt{N}})$ decay of the KRR error.

Hastie et al. (2019) use random matrix theory tools to compute the asymptotic risk when both $P,N\to\infty$ with $\frac{P}{N}\to\gamma>0$. When the training data is sampled i.i.d. from a Gaussian distribution, the variance is shown to explode at $\gamma=1$. In the same linear regression setup, Bartlett et al. (2019) establish general upper and lower bounds on the excess risk. Mei & Montanari (2019) prove that the double-descent (DD) curve also arises for random ReLU features, and adding a ridge suppresses the explosion around $\gamma=1$.

Double-descent and the effect of regularization. For the cross-entropy loss, Neyshabur et al. (2014) observed that for two-layer neural networks the test error exhibits the double-descent (DD) curve as the network width increases (without regularizers, without early stopping). For MSE and hinge losses, the DD curve was observed also in multilayer networks on the MNIST dataset (Advani & Saxe, 2017; Spigler et al., 2018). Neal et al. (2018) study the variance due to stochastic training in neural networks and find that it increases until a certain width, but then decreases down to . Nakkiran et al. (2019) establish the DD phenomenon across various models including convolutional and recurrent networks on more complex datasets (e.g. CIFAR-10, CIFAR-100).

Belkin et al. (2018, 2019) find that the DD curve is not peculiar to neural networks and observe the same for random Fourier features and decision trees. In Geiger et al. (2019), the DD curve for neural networks is related to the variance associated with the random initialization of the Neural Tangent Kernel (Jacot et al., 2018); as a result, ensembling is shown to suppress the DD phenomenon in this case, and the test error stays constant in the overparameterized regime. Recent theoretical work (d’Ascoli et al., 2020) study the same setting and derive formulas for the asymptotic error, relying on the so-called replica method.

3 Outline

The rest of this paper is organized as follows:

In Section 2, the setup (linear regression, Gaussian RF model, $\lambda$-RF predictor, and $\lambda$-KRR predictor) is introduced.

In Section 3, preliminary results on the distribution of the $\lambda$-RF model are provided: the RF predictors are Gaussian mixtures (Proposition 3.1) and the $\lambda\searrow 0$-RF model is unbiased in the overparameterized regime (Corollary 3.2). Graphical illustrations of the RF predictors in various regimes are presented (Figure 1).

In Section 6, we summarize our results and discuss potential implications and extensions.

Setup

Linear regression is a parametric model consisting of linear combinations

The data matrix $F$ is defined as the $N\times P$ matrix with entries $F_{ij}={1\over\sqrt{P}}\phi^{(j)}(x_{i})$. The minimization of (1) can be rewritten in terms of $F$ as

The optimal solution $\hat{\theta}$ is then given by

and the optimal predictor $\hat{f}=f_{\hat{\theta}}$ by

The $\lambda$-RF can be viewed as an approximation of kernel ridge predictors: observing from (4) that $\hat{f}^{(RF)}_{\lambda}$ only depends on the scalar product $K_{P}(x,x^{\prime})=\frac{1}{P}\sum_{j=1}^{P}f^{(j)}(x)f^{(j)}(x^{\prime})$ between datapoints, we see that as $P\to\infty$, $K_{P}\to K$ and hence $\hat{f}^{(RF)}_{\lambda}$ converges (Rahimi & Recht, 2008) to a kernel predictor with ridge $\lambda$ (Schölkopf et al., 1998), which we call $\lambda$-KRR predictor.

Let $\pi$ denote the joint distribution of the i.i.d. sample $f^{(1)},...,f^{(P)}$ from the centered Gaussian process with covariance kernel $K$. The risk of $\hat{f}^{(RF)}_{\lambda}$ can be decomposed into a bias-variance form as

This decomposition into the risk of the average RF predictor and of the $\mathcal{D}$-expectation of its variance will play a crucial role in the next sections. This is in contrast with the classical bias-variance decomposition in Geman et al. (1992)

where $\mathcal{D}^{\otimes N}$ denotes the joint distribution on $x_{1},...,x_{N}$, sampled i.i.d. from $\mathcal{D}$. Note that in our decomposition no probabilistic assumption is made on the data, which is fixed.

2 Additional Notation

We will denote by $\gamma=\frac{P}{N}$ the parameter-to-datapoint ratio: the underparameterized regime corresponds to $\gamma<1$, while the overparameterized regime corresponds to $\gamma\geq 1$. In order to stress the dependence on the ratio parameter $\gamma$, we write $\hat{f}^{(RF)}_{\lambda,\gamma}$ instead of $\hat{f}^{(RF)}_{\lambda}$.

First Observations

The distribution of the RF predictor features a variety of behaviors depending on $\gamma$ and $\lambda$, as displayed in Figure 1. In the underparameterized regime $P<N$, sample RF predictors induce some implicit regularization and do not interpolate the dataset (1a); at the interpolation threshold $P=N$, RF predictors interpolate the dataset but the variance explodes when there is no ridge (1b), however adding some ridge suppresses variance explosion (1c); in the overparameterized regime $P\geq N$ with large $P$, the variance vanishes thus the RF predictor converges to its average (1d). We will investigate the average RF predictor (solid lines) in detail in Section 4 and study its variance in Section 5.

We start by characterizing the distribution of the RF predictor as a Gaussian mixture:

Let $\hat{f}^{(RF)}_{\lambda,\gamma}(x)$ be the random features predictor as in (5) and let $\hat{y}=F\hat{\theta}$ be the prediction vector on training data, i.e. $\hat{y}_{i}=\hat{f}^{(RF)}_{\lambda,\gamma}(x_{i})$. The process $\hat{f}^{(RF)}_{\lambda,\gamma}$ is a mixture of Gaussians: conditioned on $F$, we have that $\hat{f}^{(RF)}_{\lambda,\gamma}$ is a Gaussian process. The mean and covariance of $\hat{f}^{(RF)}_{\lambda,\gamma}$ conditioned on $F$ are given by

The proof of Proposition 3.1 relies on the fact that $f^{(j)}$ conditioned on $\left(f^{(j)}(x_{i})\right)_{i=1,\ldots,N}$ is a Gaussian Process.

Note that (6) and (7) depend on $\lambda$ and $P$ through $\hat{y}$ and $\|\hat{\theta}\|^{2}$; in fact, as the proof shows, these identities extend to the ridgeless case $\lambda\searrow 0$. For the ridgeless case, when one is in the overparameterized regime ($P\geq N$), one can (with probability one) fit the labels $y$ and hence $\hat{y}=y$:

When $P\geq N$, the average ridgeless RF predictor is equivalent to the ridgeless KRR predictor

This corollary shows that in the overparameterized case, the ridgeless RF predictor is an unbiased estimator of the ridgeless kernel predictor. The difference between the expected loss of ridgeless RF predictor and that of the ridgeless KRR predictor is hence equal to the variance of the RF predictor. As will be demonstrated in this article, outside of this specific regime, a systematic bias appears, which reveals an implicit regularizing effect of random features.

Average Predictor

(Sketch; see Supp. Mat. for details) Set $A_{\lambda}=F(F^{T}F+\lambda I_{P})^{-1}F^{T}$. The vector of the predictions on the training set is given by $\hat{y}=A_{\lambda}y$ and the expected predictor is given by

with $W_{i}$ denoting the $i$-th column of $W=\sqrt{P}F^{T}K^{-\frac{1}{2}}$ and $F_{(i)}$ being obtained by removing the $i$-th row of $F$. The $g_{i}$’s are all within $\mathcal{O}(1/\sqrt{P})$ distance to the Stieltjes transform

(The detailed proof in the Supp. Mat. uses non-asymptotic variants of arguments found in (Bai & Wang, 2008); the constants in the $\mathcal{O}$ bounds are in particular made explicit).

As a consequence, from the above results, we obtain

Note that asymptotic forms of equations similar to the ones in the above proof appear in different settings (Dobriban & Wager, 2018; Mei & Montanari, 2019; Liu & Dobriban, 2020), related to the study of the Stieltjes transform of the product of asymptotically free random matrices.

The boundedness of $T$ is guaranteed for kernels that are translation-invariant, i.e. of the form $K(x,y)=k(\left\|x-y\right\|)$: in this case, one has $T=k(0)$.

For $\|y\|_{K^{-1}}$, under the assumption that the labels are of the form $y_{i}=f^{*}(x_{i})$ for a true regression function $f^{*}$ lying in Reproducing Kernel Hilbert Space (RKHS) $\mathcal{H}$ of the kernel $K$ (Schölkopf et al., 1998), we have $\|y\|_{K^{-1}}\leq\|f^{*}\|_{\mathcal{H}}$.

2 Risk of the Average Predictor

Variance

The following theorem allows us to bound both terms:

There are constants $c_{1},c_{2}>0$ depending on $\lambda,\gamma,T$ only such that

where $c_{3}>0$ depends on $\lambda,\gamma,T$.

2 Double Descent Curve

We now investigate the neighborhood of the frontier $\gamma=1$ between the under- and overparameterized regimes, known empirically to exhibit a double descent curve, where the test error explodes at $\gamma=1$ (i.e. when $P\approx N$) as exhibited in Figure 3.

Thanks to Theorem C.3.3, we get a lower bound on the variance of $\hat{f}^{(RF)}_{\lambda,\gamma}$:

Conclusion

Both theorems are proven using tools from random matrix theory, in particular finite-size results on the concentration of the Stieltjes transform of general Wishart matrix models. While our current proofs require the assumption that the RF model is Gaussian, it seems natural to postulate that the results and the proofs extend to more general setups, along the lines of (Louart et al., 2017; Benigni & Péché, 2019).

Finally, we investigate the ridgeless limit case $\lambda\searrow 0$. In this case, we see a sharp transition at $\gamma=1$: in the overparameterized regime $\gamma>1$, the effective ridge goes to zero, while in the underparameterized regime $\gamma<1$, it converges to a positive value. At the interpolation threshold $\gamma=1$, the variance of the $\lambda$-RF explodes, leading to the double descent curve emphasized in (Advani & Saxe, 2017; Spigler et al., 2018; Belkin et al., 2018; Nakkiran et al., 2019). We investigate this numerically and prove a lower bound yielding a plausible explanation for this phenomenon.

Thanks and Acknowledgements

The authors would like to thank Andrea Montanari, Song Mei, Lénaïc Chizat and Alessandro Rudi for the helpful discussions. Clément Hongler acknowledges support from the ERC SG CONSTAMIS grant, the NCCR SwissMAP grant, the Minerva Foundation, the Blavatnik Family Foundation, and the Latsis foundation.

References

Appendix A Experimental Details

Using the procedure above, we performed the following experiments:

A.2 MNIST experiments

A.3 Random Fourier Features

Appendix B Additional Experiments

We present the following complementary simulations:

In Section B.1, we present the distribution of the $\lambda$-RF predictor for the selected $P$ and $\lambda$.

In Section B.4, we present numerical experiments on MNIST using random Fourier features.

B.4 Average Fourier Features Predictor

The difference of the test errors of the two predictors decreases as $\gamma$ increases.

For $N=1000$, strong agreement between the two test errors is observed already for $\gamma>0.1$. We also observe that Gaussian features achieve lower (or equal) test error than the Fourier features for all $\gamma$ in our experiments.

Appendix C Proofs

Let $\hat{f}^{(RF)}_{\lambda}$ be the $\lambda$-RF predictor and let $\hat{y}=F\hat{\theta}$ be the prediction vector on training data, i.e. $\hat{y}_{i}=\hat{f}^{(RF)}_{\lambda}(x_{i})$. The process $\hat{f}^{(RF)}_{\lambda}$ is a mixture of Gaussians: conditioned on $F$, we have that $\hat{f}^{(RF)}_{\lambda}$ is a Gaussian process. The mean and covariance of $\hat{f}^{(RF)}_{\lambda}$ conditioned on $F$ are given by

Let $F=(\frac{1}{\sqrt{P}}f^{(j)}(x_{i}))_{i,j}$ be the $N\times P$ matrix of values of the random features on the training set. By definition, $\hat{f}_{\lambda}^{(RF)}=\frac{1}{\sqrt{P}}\sum_{p=1}^{P}\hat{\theta}_{p}f^{(p)}$. Conditioned on the matrix $F$, the optimal parameters $(\hat{\theta}_{p})_{p}$ are not random and $(f^{(p)})_{p}$ is still Gaussian, hence, conditioned on the matrix $F$, the process $\hat{f}_{\lambda}^{(RF)}$ is a mixture of Gaussians. Moreover, conditioned on the matrix $F$, for any $p,p^{\prime}$, $f^{(p)}$ and $f^{(p^{\prime})}$ remain independent, hence

C.2 Generalized Wishart Matrix

Setup. In this section, we consider a fixed deterministic matrix $K$ of size $N\times N$ which is diagonal positive semi-definite, with eigenvalues $d_{1},\ldots,d_{N}$. We also consider a $P\times N$ random matrix $W$ with i.i.d. standard Gaussian entries.

where $K$ is a fixed positive semi-definite matrix.

Our first lemma implies that the Stieljes transform concentrates around its mean as $N$ and $P$ go to infinity with $\gamma=\frac{P}{N}$ fixed.

where $\mathbf{c}$ depends on $z$, $\gamma$, and $m$ only.

is a rank one perturbation of the matrix $B_{(i)}(z)$, by the Sherman–Morrison’s formula, the inverse of $B(z)$ is given by:

where we used the cyclic property of the trace. We can now bound this difference:

where $\nu_{j}$ are the eigenvalues of $\frac{1}{P}W_{(i)}K_{(i)}W_{(i)}^{T}$.

is a martingale difference sequence. Hence, by Burkholder’s inequality, there exists a positive constant $K_{m}$ such that

The following lemma, which is reminiscent of Lemma 4.5 in (Au et al., 2018), is a consequence of Wick’s formula for Gaussian random variables and is key to prove Lemma C.4.

If $A^{(1)},\ldots,A^{(k)}$ are $k$ square random matrices of size $P$ independent from a standard Gaussian vector $w$ of size $P$,

where $:\boldsymbol{P}_{2}(2k):$ is the subset of partitions $p$ in $\boldsymbol{P}_{2}(2k)$ for which $\left\{2j-1,2j\right\}$ is not a block of $p$ for any $j\in\{1,\ldots,k\}$.

Expanding the left-hand side of Equation (12), we obtain:

hence, interchanging the order of summation, we recover the left-hand side of Equation (12):

We now prove Equation (13). Expanding the product, the left-hand side is equal to:

Expanding the product and the trace, and using Wick’s equation, we obtain: a

where $p_{I}$ is the partition composed of blocks of size $2$ given by $\{2l,2l+1\}$ with $l\notin I$ and the rest of the indices contained in a single block. Interchanging the order of summation, we get:

Since $\left[\sum_{I\subset\{1,\ldots,k\},\!p\leq p_{I}}(-1)^{\#I}\right]=\delta_{\{I\subset[k],p\leq p_{I}\}=\{\{1,\ldots,k\}\}}$ and $\{I\subset[k],p\leq p_{I}\}=\{\{1,\ldots,k\}\}$ if and only if $p\in:\!\!\!\boldsymbol{P}_{2}(2k)\!\!\!:$, interchanging a last time the order of summation, we recover the left-hand side of Equation (13):

The random function $g_{i}(z)$ satisfies:

where $\mathbf{c_{0}}$, $\mathbf{c_{1}}$, $\mathbf{c_{2}}$, and $\mathbf{c_{3}}$ depend on $\gamma$ and $z$ only.

Thus, the variance of $g_{i}(z)$ is given by:

We bound the second term using the concentration of the Stieljes transform (Lemma C.2): it is bounded by $\frac{8\mathbf{c}}{P^{2}}$. The first term is bounded using the second assertion of Lemma C.3. Using the symmetry of $B_{(i)}(z)$, the partitions in $:\boldsymbol{P}_{2}(4):$ yield two different terms, namely:

In the next proposition we show that the Stieltjes transform $m_{P}(z)$ is close in expectation to the solution of a fixed point equation.

From the proof of Lemma C.2, recall that $B^{-1}(z)=B_{(i)}^{-1}(z)-\frac{d_{i}}{P}\frac{1}{1+\frac{d_{i}}{P}w_{i}^{T}B_{(i)}^{-1}(z)w_{i}}B_{(i)}^{-1}(z)w_{i}w_{i}^{T}B_{(i)}^{-1}(z),$ hence:

Thus, the Stieljes transform $m_{P}(z)$ satisfies the following equation $P=\sum_{i=1}^{N}\frac{d_{i}g_{i}(z)}{1+d_{i}g_{i}(z)}-zPm_{P}(z),$ or equivalently

To prove that $0\in\psi(\mathcal{C}_{z})$ we proceed with a geometrical reasoning: the image $\psi(\mathcal{C}_{z})$ is (one of) the region of the plane confined by $\psi\left(\partial\mathcal{C}_{z}\right)$, so we only need to “draw” $\psi\left(\partial\mathcal{C}_{z}\right)$ and show that belongs to the “good” connected component confined by it.

We observe that, for every $t\in\mathcal{C}_{z}$, the derivative of $f$ has negative real part:

where we concluded the last inequality by using that $\Re(z)\leq 0$, $\Re(t)\geq 0$, $\Im(z)\Im(t)\geq 0$ and $\Re(zt^{2})\leq 0$. Thus, since for no point $t\in\mathcal{C}_{z}$ has $f_{z}^{\prime}(t)=1$, any fixed point of $f_{z}$ is a simple fixed point.

We now proceed to show the uniqueness of the fixed point in the region $\mathcal{C}_{z}$. Suppose there are two fixed points $t_{1}$ and $t_{2}$, then

We provide a lower bound on the norm of the fixed point:

C.3 Ridge

Note that the matrix $A_{\lambda}$ defined in the proof sketch of Theorem 4.1 in the main text is given by $A_{\lambda}=A(-\lambda)$.

Since the distribution of $W$ is invariant under orthogonal transformations, by applying a change of basis, in order to prove Inequality (17), we may assume that $K$ is diagonal with diagonal entries $d_{1},\ldots,d_{N}$. Denoting $w_{1},\ldots,w_{N}$ the columns of $W$, for any $i,j=1,\ldots,N$,

Consider a diagonal term $(A(z))_{ii}$. From Equation (15), we get

For the second term, using the same arguments as for the proof of Proposition C.5, we have:

For a general vector $v$, the $K^{-1}$-norm $\left\|v\right\|_{K^{-1}}$ is equal to the norm mininum Hilbert norm (for the RKHS associated to the kernel $K$) interpolating function:

Indeed the minimal interpolating function is the kernel regression given by $f^{(K)}(\cdot)=K(\cdot,X)K(X,X)^{-1}v$ which has norm (writing $\beta=K^{-1}v$):

We can now bound the two norms $\left\|v\right\|_{K^{-1}}$ and $\left\|w\right\|_{K^{-1}}$. For $v=K(x,X)$, we have

since $K(x,\cdot)$ is an interpolating function for $v$.

Hence, if $f^{*}$ is the true function, by the triangular inequality,

C.3.2 Properties of the effective ridge

Differentiating both sides of Equation (20),

C.3.3 Variance of the predictor

This section is dedicated to the proof of the variance bound of Theorem 5.1 of the paper:

Theorem 5.1 There are constants $c_{1},c_{2}>0$ depending on $\lambda,\gamma,T$ only such that

where $c_{3}>0$ depends on $\lambda,\gamma,T$.

Let us now consider the second term in the r.h.s. of (23) . Using the fact that $\|B_{(i)}\|_{op}\geq\frac{1}{\lambda}$, we get

where we have used the fact that the second moment of a $\chi^{2}(N)$ distribution is $N(N+2).$ Together, we obtain

There exists a constant $c_{1}>0$ (depending on $\lambda,\gamma,T$ only) such that the variance of the estimator is bounded by

As in the proof of Theorem C.8, with the right change of basis, we may assume the Gram matrix $K(X,X)$ to be diagonal.

We first express the covariances of $\hat{y}=A(-\lambda)y$. Using Proposition Proposition C.12, for $i\neq j$ we have

and let $D$ the diagonal matrix with entries

which implies that $\|D\|_{op}\leq\frac{c^{\prime}_{1}\|y\|_{K^{-1}}^{2}}{P}$. As a result

where we used Inequality (21). This yields the result with $c_{1}=3c^{\prime}_{1}$.∎

For $\gamma,\lambda>0$ there exists a constant $c_{2}>0$ depending on $\lambda,\gamma,T$ only such that

As in the proof of Theorem C.8, with the right change of basis, we may assume the Gram matrix $K(X,X)$ to be diagonal. Recall that $\hat{\theta}=\frac{1}{\sqrt{P}}\left(\frac{1}{P}WK(X,X)W^{T}+\lambda I_{N}\right){}^{-1}WK(X,X)^{\frac{1}{2}}y$, thus we have:

where $A^{\prime}(-\lambda)$ is the derivative of

with respect to $z$ evaluated at $-\lambda$. Let

where $c_{3}>0$ depends on $\lambda,\gamma,T$.

Using the bias/variance decomposition, Corollary C.9, and the bound on the variance of the predictor, we obtain

C.3.5 Double descent curve