A random matrix analysis of random Fourier features: beyond the Gaussian kernel, a precise phase transition, and the corresponding double descent^*

We consider the RFF model in the more realistic large n, p, N limit. While, in this setting, the RFF empirical Gram matrix does not converge to the Gaussian kernel matrix, we can characterize its behavior as n, p, N → ∞ and provide asymptotic performance guarantees for RFF on large-scale problems. We also identify a phase transition as a function of the ratio N/n, including the corresponding double descent phenomenon. In more detail, our contributions are the following.

• (a)  
  We provide a precise characterization of the asymptotics of the RFF empirical Gram matrix, in the large n, p, N limit (theorem 1). This is accomplished by constructing a deterministic equivalent for the resolvent of the RFF Gram matrix. Based on this, the behavior of the RFF model is (asymptotically) accessible through a fixed-point equation, that can be interpreted in terms of an angle-like correction induced by the non-trivial large n, p, N limit (relative to the N → ∞ alone limit).
• (b)  
  We derive the asymptotic training and test mean squared errors (MSEs) of RFF ridge regression, as a function of the ratio N/n, regularization penalty λ, training as well as test sets (theorems 2 and 3, respectively). We identify precisely the under-to over-parameterization phase transition, as a function of the relative model complexity N/n; we prove the existence of a 'singular' peak of test error at the N/n = 1/2 boundary; and we characterize the corresponding double descent behavior. Importantly, our results are valid with almost no specific assumption on the data distribution. This is a significant improvement over existing double descent analyses, which fundamentally rely on the knowledge of the data distribution (often assumed to be multivariate Gaussian for simplicity) [12, 18].
• (c)  
  We provide a detailed empirical evaluation of our theoretical results, demonstrating that the theory closely matches empirical results on a range of real-world data sets (sections 3 and 4). This includes the correction due to the large n, p, N setting, sharp transitions (as a function of N/n) in the aforementioned angle-like quantities, and the corresponding double descent test curves. This also includes an evaluation of the impact of training-test similarity and the effect of different data sets, thus confirming, as stated in (ii), that (unlike in prior work) the phase transition and double descent curve hold much more generally with respect to the data distribution.

Here, we provide a brief review of related previous efforts.

Random features and limiting kernels. In most RFF work [19–22], non-asymptotic bounds are given, on the number of random features N needed for a predefined approximation error, for a given kernel matrix with fixed n, p. A more recent line of work [1, 23–25] has focused on the over-parameterized N → ∞ limit of large neural networks by studying the corresponding NTKs. Here, we position ourselves in the more practical regime where n, p, N are all large and comparable, and we provide asymptotic performance guarantees that better fit large-scale problems compared to the large-N-alone analysis.

Random matrix theory. From a random matrix theory perspective, nonlinear Gram matrices of the type ${\mathbf{\Sigma }}_{\mathbf{X}}^{\mathsf{\text{T}}}{\mathbf{\Sigma }}_{\mathbf{X}}$ have recently received an unprecedented research interests, due to their close connection to neural networks [26–29], with a particular focus on the associated eigenvalue distribution. Here we propose a deterministic equivalent [30, 31] analysis for the resolvent matrix that provides access, not only to the eigenvalue distribution, but also to the regression error of central interest in this article. While most existing deterministic equivalent analyses are performed on linear models, here we focus on the nonlinear RFF model. From a technical perspective, the most relevant work is [12, 15]. We improve their results by considering generic data model on the popular RFF model.

Statistical mechanics of learning. A long history of connections between statistical mechanics and machine learning models (such as neural networks) exists, including a range of techniques to establish generalization bounds [2–5], and recently there has been renewed interest [7, 8, 32–34]. Their relevance to our results lies in the use of the so-called thermodynamic limit (akin to the large n, p, N limit), rather than the classical limits more commonly used in statistical learning theory, in which case uniform convergence bounds and related techniques can be applied.

Double descent in large-scale learning systems. The large n, N asymptotics of statistical models has received considerable research interests in the machine learning community [18, 35], resulting in a (somehow) counterintuitive phenomenon referred to as the 'double descent'. Instead of focusing on different 'phases of learning' [2–5, 7], the 'double descent' phenomenon focuses on an empirical manifestation of the phase boundary and refers to the empirical observations of the test error curve as a function of the model complexity, which differs from the usual textbook description of the bias-variance tradeoff [10, 11, 36, 37]. Theoretical investigation into this phenomenon mainly focuses on various regression models [12, 18, 38–41]. In most cases, quite specific (and rather strong) assumptions are imposed on the input data distribution. In this respect, our work extends the analysis in [12] to handle the RFF model and its phase structure on real-world data sets.

Throughout this article, we follow the convention of denoting scalars by lowercase, vectors by lowercase boldface, and matrices by uppercase boldface letters. In addition, the notation (⋅) ^T denotes the transpose operator; the norm ||⋅|| is the Euclidean norm for vectors and the spectral or operator norm for matrices; and $\;\stackrel{a.s.}{\to }\;$ stands for almost sure convergence of random variables.

Our main results on the asymptotic behavior of the RFF resolvent matrix, as well as of the training MSE and testing MSE of RFF ridge regression are presented in section 2, with detailed proofs deferred to the appendix. In section 3, we provide a detailed empirical evaluation of our main results; and in section 4, we provide additional empirical evaluation on real-world data, illustrating the practical effectiveness of the proposed analysis. Concluding remarks are placed in section 5.

To start, we observe that the training MSE, E_train, in (5), can be written as

$$\begin{equation}{E}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}=\frac{{\lambda }^{2}}{n}{\Vert}\mathbf{Q}(\lambda )\mathbf{y}{{\Vert}}^{2}=-\frac{{\lambda }^{2}}{n}{\mathbf{y}}^{\mathsf{\text{T}}}\partial \mathbf{Q}(\lambda )\mathbf{y}/\partial \lambda ,\end{equation} \tag{ 6 }$$

which depends on the quadratic form y^T Q(λ)y of

$$\begin{equation}\mathbf{Q}(\lambda )\equiv {\left(\frac{1}{n}{\mathbf{\Sigma }}_{\mathbf{X}}^{\mathsf{\text{T}}}{\mathbf{\Sigma }}_{\mathbf{X}}+\lambda {\mathbf{I}}_{n}\right)}^{-1}\in {\mathbb{R}}^{n\times n},\end{equation} \tag{ 7 }$$

the so-called resolvent of $\frac{1}{n}{\mathbf{\Sigma }}_{\mathbf{X}}^{\mathsf{\text{T}}}{\mathbf{\Sigma }}_{\mathbf{X}}$ (also denoted Q when there is no ambiguity) with λ > 0. To see this, from (5) we have ${E}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}=\frac{1}{n}{\Vert}\mathbf{y}-\frac{1}{n}{\mathbf{\Sigma }}_{\mathbf{X}}^{\mathsf{\text{T}}}{\mathbf{\Sigma }}_{\mathbf{X}}{(\frac{1}{n}{\mathbf{\Sigma }}_{\mathbf{X}}^{\mathsf{\text{T}}}{\mathbf{\Sigma }}_{\mathbf{X}}+\lambda {\mathbf{I}}_{n})}^{-1}\mathbf{y}{{\Vert}}^{2}=\frac{{\lambda }^{2}}{n}{\Vert}\mathbf{Q}(\lambda )\mathbf{y}{{\Vert}}^{2}=-\frac{{\lambda }^{2}}{n}{\mathbf{y}}^{\mathsf{\text{T}}}\frac{\partial \mathbf{Q}(\lambda )}{\partial \lambda }\mathbf{y}$, with $\frac{\partial \mathbf{Q}(\lambda )}{\partial \lambda }=-{\mathbf{Q}}^{2}(\lambda )$.

In order to assess the asymptotic training MSE, it thus suffices to find a deterministic equivalent for Q(λ), that is, a deterministic matrix that captures the asymptotic behavior of the latter. One possibility is the expectation ${\mathbb{E}}_{\mathbf{W}}[\mathbf{Q}(\lambda )]$. Informally, if the training MSE E_train (that is random due to random W for given X, y) is 'close to' some deterministic quantity ${\bar{E}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}$, in the large n, p, N limit, then ${\bar{E}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}$ must have the same limit as ${\mathbb{E}}_{\mathbf{W}}[{E}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}]=-\frac{{\lambda }^{2}}{n}\partial {\mathbf{y}}^{\mathsf{\text{T}}}{\mathbb{E}}_{\mathbf{W}}[\mathbf{Q}(\lambda )]\mathbf{y}/\partial \lambda$ for n, p, N → ∞. However, ${\mathbb{E}}_{\mathbf{W}}[\mathbf{Q}]$ involves integration (with no closed-form due to the matrix inverse), and it is not a convenient quantity with which to work. Our objective is to find an asymptotic 'alternative' for ${\mathbb{E}}_{\mathbf{W}}[\mathbf{Q}]$ that is (i) close to ${\mathbb{E}}_{\mathbf{W}}[\mathbf{Q}]$ in the large n, p, N → ∞ limit and (ii) numerically more accessible.

In the following theorem, we introduce an asymptotic equivalent for ${\mathbb{E}}_{\mathbf{W}}[\mathbf{Q}]$. Instead of being directly related to the Gaussian kernel K_X = K_cos + K_sin as suggested by (2) in the large-N-only limit, it depends on the two components K_cos, K_sin in a more involved manner. Importantly, the proposed equivalent $\bar{\mathbf{Q}}$ can be numerically evaluated by running simple fixed-point iterations involving K_cos and K_sin.

Theorem 1. (Asymptotic equivalent for ${\mathbb{E}}_{\mathbf{W}}$[Q]). Under assumption 1, for Q defined in (7) and λ > 0, we have, as n → ∞

$$\begin{equation*}{\Vert}{\mathbb{E}}_{\mathbf{W}}[\mathbf{Q}]-\bar{\mathbf{Q}}{\Vert}\to 0\end{equation*}$$

for $\bar{\mathbf{Q}}\equiv {\left(\frac{N}{n}(\frac{{\mathbf{K}}_{\mathrm{cos}}}{1+{\delta }_{\mathrm{cos}}}+\frac{{\mathbf{K}}_{\mathrm{sin}}}{1+{\delta }_{\mathrm{sin}}})+\lambda {\mathbf{I}}_{n}\right)}^{-1}$, ${\mathbf{K}}_{\mathrm{cos}}\equiv {\mathbf{K}}_{\mathrm{cos}}(\mathbf{X},\mathbf{X}),{\mathbf{K}}_{\mathrm{sin}}\equiv {\mathbf{K}}_{\mathrm{sin}}(\mathbf{X},\mathbf{X})\in {\mathbb{R}}^{n\times n}$ and

$$\begin{align}\hfill \begin{aligned}\hfill {\mathbf{K}}_{\mathrm{cos}}{(\mathbf{X},{\mathbf{X}}^{\prime })}_{ij}& ={\text{e}}^{-\frac{{\Vert}{\mathbf{x}}_{i}{{\Vert}}^{2}+{\Vert}{\mathbf{x}}_{j}^{\prime }{{\Vert}}^{2}}{2}}\enspace \mathrm{cosh}({\mathbf{x}}_{i}^{\mathsf{\text{T}}}{\mathbf{x}}_{j}^{\prime }),\hfill \\ \hfill {\mathbf{K}}_{\mathrm{sin}}{(\mathbf{X},{\mathbf{X}}^{\prime })}_{ij}& ={\text{e}}^{-\frac{{\Vert}{\mathbf{x}}_{i}{{\Vert}}^{2}+{\Vert}{\mathbf{x}}_{j}^{\prime }{{\Vert}}^{2}}{2}}\enspace \mathrm{sinh}({\mathbf{x}}_{i}^{\mathsf{\text{T}}}{\mathbf{x}}_{j}^{\prime }),\hfill \end{aligned}\end{align} \tag{ 8 }$$

where (δ_cos, δ_sin) is the unique positive solution to

$$\begin{equation}{\delta }_{\mathrm{cos}}=\frac{1}{n}\enspace \mathrm{t}\mathrm{r}({\mathbf{K}}_{\mathrm{cos}}\bar{\mathbf{Q}}),\quad {\delta }_{\mathrm{sin}}=\frac{1}{n}\enspace \mathrm{t}\mathrm{r}({\mathbf{K}}_{\mathrm{sin}}\bar{\mathbf{Q}}).\end{equation} \tag{ 9 }$$

Proof. See appendix A. □

Remark 1. (Lower and upper bounds). Since

$$\begin{equation}\frac{{\mathbf{K}}_{\mathbf{X}}}{1+\mathrm{max}({\delta }_{\mathrm{cos}},{\delta }_{\mathrm{sin}})}{\preceq}\frac{{\mathbf{K}}_{\mathrm{cos}}}{1+{\delta }_{\mathrm{cos}}}+\frac{{\mathbf{K}}_{\mathrm{sin}}}{1+{\delta }_{\mathrm{sin}}}{\preceq}\frac{{\mathbf{K}}_{\mathbf{X}}}{1+\mathrm{min}({\delta }_{\mathrm{cos}},{\delta }_{\mathrm{sin}})}\end{equation} \tag{ 10 }$$

in the positive definite order, for K_X ≡ K_cos + K_sin the Gaussian kernel, $\frac{{\mathbf{K}}_{\mathrm{cos}}}{1+{\delta }_{\mathrm{cos}}}+\frac{{\mathbf{K}}_{\mathrm{sin}}}{1+{\delta }_{\mathrm{sin}}}$ is therefore positive definite, if x₁ ..., x_n are all distinct; see theorem 2.18 in [42].

Remark 2. (Correction to large-N behavior). Taking N/n → ∞, one has δ_cos → 0, δ_sin → 0 so that

$$\begin{equation}\frac{{\mathbf{K}}_{\mathrm{cos}}}{1+{\delta }_{\mathrm{cos}}}+\frac{{\mathbf{K}}_{\mathrm{sin}}}{1+{\delta }_{\mathrm{sin}}}\to {\mathbf{K}}_{\mathrm{cos}}+{\mathbf{K}}_{\mathrm{sin}}={\mathbf{K}}_{\mathbf{X}}\quad \text{and}\quad \bar{\mathbf{Q}}\to {\left(\frac{N}{n}{\mathbf{K}}_{\mathbf{X}}+\lambda {\mathbf{I}}_{n}\right)}^{-1}\sim \frac{n}{N}{\mathbf{K}}_{\mathbf{X}}^{-1},\end{equation} \tag{ 11 }$$

for λ > 0 independent of N, n, in accordance with the classical large-N-only prediction. In this sense, the pair (δ_cos, δ_sin) introduced in theorem 1 accounts for the 'correction' due to the non-trivial n/N, as opposed to the N → ∞ alone analysis. Also, when the number of features N is large (i.e. as N/n → ∞), the regularization effect of λ flattens out and $\bar{\mathbf{Q}}$ behaves like (a scaled version of) the inverse Gaussian kernel matrix ${\mathbf{K}}_{\mathbf{X}}^{-1}$ (that is well-defined for distinct x₁ ..., x_n ).

Remark 3. (Geometric interpretation). Since $\bar{\mathbf{Q}}$ shares the same eigenspace with $\frac{{\mathbf{K}}_{\mathrm{cos}}}{1+{\delta }_{\mathrm{cos}}}+\frac{{\mathbf{K}}_{\mathrm{sin}}}{1+{\delta }_{\mathrm{sin}}}$, one can geometrically interpret (δ_cos, δ_sin) as a sort of 'angle' between the eigenspaces of K_cos, K_sin and that of $\frac{{\mathbf{K}}_{\mathrm{cos}}}{1+{\delta }_{\mathrm{cos}}}+\frac{{\mathbf{K}}_{\mathrm{sin}}}{1+{\delta }_{\mathrm{sin}}}$. For fixed n, as N → ∞, one has $\frac{1}{N}{\sum }_{t=1}^{N}\mathrm{cos}({\mathbf{X}}^{\mathsf{\text{T}}}{\mathbf{w}}_{t})\mathrm{cos}({\mathbf{w}}_{t}^{\mathsf{\text{T}}}\mathbf{X})\to {\mathbf{K}}_{\mathrm{cos}}$, $\frac{1}{N}{\sum }_{t=1}^{N}\mathrm{sin}({\mathbf{X}}^{\mathsf{\text{T}}}{\mathbf{w}}_{t})\mathrm{sin}({\mathbf{w}}_{t}^{\mathsf{\text{T}}}\mathbf{X})\to {\mathbf{K}}_{\mathrm{sin}}$, the eigenspaces of which are 'orthogonal' to each other, so that δ_cos, δ_sin → 0. On the other hand, as N, n → ∞, the eigenspaces of K_cos and K_sin 'intersect' with each other, captured by the non-trivial (δ_cos, δ_sin).

Theorem 1 provides an asymptotically more tractable approximation of ${\mathbb{E}}_{\mathbf{W}}[\mathbf{Q}]$. Together with some additional concentration arguments (e.g. from theorem 2 in [15]). this permits us to provide a complete description of the limiting behavior of the random bilinear form a ^T Qb, for $\mathbf{a},\mathbf{b}\in {\mathbb{R}}^{n}$ of bounded Euclidean norms, in such a way that ${\mathbf{a}}^{\mathsf{\text{T}}}\mathbf{Q}\mathbf{b}-{\mathbf{a}}^{\mathsf{\text{T}}}\bar{\mathbf{Q}}\mathbf{b}\;\stackrel{a.s.}{\to }\;0$, as n, p, N → ∞. This, together with the fact that ${E}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}=\frac{{\lambda }^{2}}{n}{\mathbf{y}}^{\mathsf{\text{T}}}\mathbf{Q}{(\lambda )}^{2}\mathbf{y}=-\frac{{\lambda }^{2}}{n}{\mathbf{y}}^{\mathsf{\text{T}}}\partial \mathbf{Q}(\lambda )\mathbf{y}/\partial \lambda$, leads to the following result on the asymptotic training error.

Theorem 2. (Asymptotic training performance). Under assumption 1, for a given training set (X, y) and training MSE, E_train defined in (5), as n → ∞

$$\begin{align*}\hfill {E}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}-{\bar{E}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}\;\stackrel{a.s.}{\to }\;0,\quad {\bar{E}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}& =\frac{{\lambda }^{2}}{n}{\Vert}\bar{\mathbf{Q}}\mathbf{y}{{\Vert}}^{2}+\frac{N}{n}\frac{{\lambda }^{2}}{{n}^{2}}\left[\begin{matrix}\hfill \frac{\mathrm{t}\mathrm{r}(\bar{\mathbf{Q}}{\mathbf{K}}_{\mathrm{cos}}\bar{\mathbf{Q}})}{{(1+{\delta }_{\mathrm{cos}})}^{2}}\hfill & \hfill \frac{\mathrm{t}\mathrm{r}(\bar{\mathbf{Q}}{\mathbf{K}}_{\mathrm{sin}}\bar{\mathbf{Q}})}{{(1+{\delta }_{\mathrm{sin}})}^{2}}\hfill \end{matrix}\right]\hfill \\ \hfill & \quad \times \mathbf{\Omega }\left[\begin{matrix}\hfill {\mathbf{y}}^{\mathsf{\text{T}}}\bar{\mathbf{Q}}{\mathbf{K}}_{\mathrm{cos}}\bar{\mathbf{Q}}\mathbf{y}\hfill \\ \hfill {\mathbf{y}}^{\mathsf{\text{T}}}\bar{\mathbf{Q}}{\mathbf{K}}_{\mathrm{sin}}\bar{\mathbf{Q}}\mathbf{y}\hfill \end{matrix}\right]\hfill \end{align*}$$

for $\bar{\mathbf{Q}}$ defined in theorem 1 and

$$\begin{equation}{\mathbf{\Omega }}^{-1}\equiv {\mathbf{I}}_{2}-\frac{N}{n}\left[\begin{matrix}\hfill \frac{1}{n}\frac{\mathrm{t}\mathrm{r}(\bar{\mathbf{Q}}{\mathbf{K}}_{\mathrm{cos}}\bar{\mathbf{Q}}{\mathbf{K}}_{\mathrm{cos}})}{{(1+{\delta }_{\mathrm{cos}})}^{2}}\hfill & \hfill \frac{1}{n}\frac{\mathrm{t}\mathrm{r}(\bar{\mathbf{Q}}{\mathbf{K}}_{\mathrm{cos}}\bar{\mathbf{Q}}{\mathbf{K}}_{\mathrm{sin}})}{{(1+{\delta }_{\mathrm{sin}})}^{2}}\hfill \\ \hfill \frac{1}{n}\frac{\mathrm{t}\mathrm{r}(\bar{\mathbf{Q}}{\mathbf{K}}_{\mathrm{cos}}\bar{\mathbf{Q}}{\mathbf{K}}_{\mathrm{sin}})}{{(1+{\delta }_{\mathrm{cos}})}^{2}}\hfill & \hfill \frac{1}{n}\frac{\mathrm{t}\mathrm{r}(\bar{\mathbf{Q}}{\mathbf{K}}_{\mathrm{sin}}\bar{\mathbf{Q}}{\mathbf{K}}_{\mathrm{sin}})}{{(1+{\delta }_{\mathrm{sin}})}^{2}}\hfill \end{matrix}\right].\end{equation} \tag{ 12 }$$

Proof. See appendix B. □

Remark 4. (First- and second-order corrections). Since ${E}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}=\frac{{\lambda }^{2}}{n}{\mathbf{y}}^{\mathsf{\text{T}}}{\mathbf{Q}}^{2}\mathbf{y}$, we can see in the expression of ${\bar{E}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}$ that there is not only a first-order (large n, p, N) correction in the first $\frac{{\lambda }^{2}}{n}{\Vert}\bar{\mathbf{Q}}\mathbf{y}{{\Vert}}^{2}$ term (which is different than $\frac{{\lambda }^{2}}{n}{\Vert}\mathbf{Q}\mathbf{y}{{\Vert}}^{2}$), but there is also a second-order correction, appearing in the form of $\bar{\mathbf{Q}}{\mathbf{K}}_{\sigma }\bar{\mathbf{Q}}$ or $\bar{\mathbf{Q}}{\mathbf{K}}_{\sigma }\bar{\mathbf{Q}}{\mathbf{K}}_{\sigma }$ for σ ∈ { cos,  sin }, as in the second term. This has a similar interpretation to remark 3, where the pair (δ_cos, δ_sin) in $\bar{\mathbf{Q}}$ is (geometrically) interpreted as the eigenspace 'intersection' due to a non-vanishing n/N. In particular, taking N/n → ∞, we have $\bar{\mathbf{Q}}\sim \frac{n}{N}{\mathbf{K}}_{\mathbf{X}}^{-1}$, Ω → I₂ so that ${\bar{E}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}=0$ and the model interpolates the entire training set, as expected.

One can show that (i) for a given n and λ > 0, ${\bar{E}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}$ decreases as the model size N increases; and (ii) for a given ratio N/n, ${\bar{E}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}$ increases as the regularization penalty λ grows large.

Theorem 2 holds without any restriction on the training set, (X, y), except for assumption 1, since only the randomness of W is involved, and thus one can simply treat (X, y) as known in this result. This is no longer the case for the test error. Intuitively, the test data $\hat{\mathbf{X}}$ cannot be chosen arbitrarily, and one must ensure that the test data 'behave' statistically like the training data, in some 'well-controlled' manner, so that the test MSE is asymptotically deterministic and bounded as $n,\hat{n},p,N\to \infty$. Following this intuition, we work under the following assumption.

Assumption 2. (Data as concentrated random vectors [43]). The training data ${\mathbf{x}}_{i}\in {\mathbb{R}}^{p},i\in \left\{1,\dots ,n\right\}$, are independently drawn (non-necessarily uniformly) from one of K > 0 distribution classes ⁴ μ₁, ..., μ_K . There exist constants C, η, q > 0 such that for any x_i ∼ μ_k , k ∈ {1, ..., K} and any one-Lipschitz function $f:{\mathbb{R}}^{p}\to \mathbb{R}$, we have

$$\begin{equation}\mathbb{P}\left(\left\vert f({\mathbf{x}}_{i})-\mathbb{E}[f({\mathbf{x}}_{i})]\right\vert \geqslant t\right)\leqslant C{e}^{-{(t/\eta )}^{q}},\quad t\geqslant 0.\end{equation} \tag{ 13 }$$

The test data ${\hat{\mathbf{x}}}_{i}\sim {\mu }_{k}$, $i\in \left\{1,\dots ,\hat{n}\right\}$ are mutually independent, but may depend on training data X and that ${\Vert}\mathbb{E}[\sigma (\mathbf{W}\mathbf{X})-\sigma (\mathbf{W}\hat{\mathbf{X}})]{\Vert}=O(\sqrt{n})$ for σ ∈ { cos,  sin }.

To facilitate the discussion of the phase transition and the double descent, we do not assume independence between training and test data (but we do assume independence between different data vectors within X and $\hat{\mathbf{X}}$). In this respect, assumption 2 is weaker than the classical i.i.d. assumption, and it permits us to illustrate the impact of training-test similarity on the model performance (section 4.2).

A first example of concentrated random vectors satisfying (13) is the multivariate Gaussian vector $\mathcal{N}(\mathbf{0},{\mathbf{I}}_{p})$ [44]. Moreover, since the concentration property in (13) is stable over Lipschitz transformations [43], it holds, for any one-Lipschitz mapping $g:{\mathbb{R}}^{d}\to {\mathbb{R}}^{p}$ and $\mathbf{z}\sim \mathcal{N}(\mathbf{0},{\mathbf{I}}_{d})$, that g(z) also satisfies (13). In this respect, assumption 2, although seemingly quite restrictive, represents a large family of 'generative models', including notably the 'fake images' generated by modern generative adversarial networks that are, by construction, Lipschitz transformations of large random Gaussian vectors [45, 46]. As such, from a practical consideration, assumption 2 provides a more realistic and flexible statistical model for real-world data.

With assumption 2, we have the following result on the asymptotic test error.

Theorem 3. (Asymptotic test performance). Under assumptions 1 and 2, we have, for test MSE E_test defined in (5) and test data $(\hat{\mathbf{X}},\hat{\mathbf{y}})$ satisfying ${\mathrm{limsup}}_{\hat{n}}{\Vert}\hat{\mathbf{X}}{\Vert}< \infty$, ${\mathrm{limsup}}_{\hat{n}}{\Vert}\hat{\mathbf{y}}{{\Vert}}_{\infty }< \infty$ with $\hat{n}/n\in (0,\infty )$ that, as n → ∞

$$\begin{align*}\hfill {E}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}-{\bar{E}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}\;\stackrel{a.s.}{\to }\;0,\quad {\bar{E}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}& =\frac{1}{\hat{n}}{\Vert}\hat{\mathbf{y}}-\frac{N}{n}\hat{\mathbf{\Phi }}\bar{\mathbf{Q}}\mathbf{y}{{\Vert}}^{2}+\frac{{N}^{2}}{{n}^{2}\hat{n}}\left[\begin{matrix}\hfill \frac{{{\Theta}}_{\mathrm{cos}}}{{(1+{\delta }_{\mathrm{cos}})}^{2}}\hfill & \hfill \frac{{{\Theta}}_{\mathrm{sin}}}{{(1+{\delta }_{\mathrm{sin}})}^{2}}\hfill \end{matrix}\right]\hfill \\ \hfill & \quad \times \mathbf{\Omega }\left[\begin{matrix}\hfill {\mathbf{y}}^{\mathsf{\text{T}}}\bar{\mathbf{Q}}{\mathbf{K}}_{\mathrm{cos}}\bar{\mathbf{Q}}\mathbf{y}\hfill \\ \hfill {\mathbf{y}}^{\mathsf{\text{T}}}\bar{\mathbf{Q}}{\mathbf{K}}_{\mathrm{sin}}\bar{\mathbf{Q}}\mathbf{y}\hfill \end{matrix}\right]\hfill \end{align*}$$

for Ω defined in (12),

$$\begin{equation}{{\Theta}}_{\sigma }=\frac{1}{N}\enspace \mathrm{t}\mathrm{r}\enspace {\mathbf{K}}_{\sigma }(\hat{\mathbf{X}},\hat{\mathbf{X}})+\frac{N}{n}\frac{1}{n}\enspace \mathrm{t}\mathrm{r}\enspace \bar{\mathbf{Q}}{\hat{\mathbf{\Phi }}}^{\mathsf{\text{T}}}\hat{\mathbf{\Phi }}\bar{\mathbf{Q}}{\mathbf{K}}_{\sigma }-\frac{2}{n}\enspace \mathrm{t}\mathrm{r}\enspace \bar{\mathbf{Q}}{\hat{\mathbf{\Phi }}}^{\mathsf{\text{T}}}{\mathbf{K}}_{\sigma }(\hat{\mathbf{X}},\mathbf{X}),\quad \sigma \in \left\{\mathrm{cos},\mathrm{sin}\right\},\end{equation} \tag{ 14 }$$

and $\mathbf{\Phi }\equiv \frac{{\mathbf{K}}_{\mathrm{cos}}}{1+{\delta }_{\mathrm{cos}}}+\frac{{\mathbf{K}}_{\mathrm{sin}}}{1+{\delta }_{\mathrm{sin}}}$, $\hat{\mathbf{\Phi }}\equiv \frac{{\mathbf{K}}_{\mathrm{cos}}(\hat{\mathbf{X}},\mathbf{X})}{1+{\delta }_{\mathrm{cos}}}+\frac{{\mathbf{K}}_{\mathrm{sin}}(\hat{\mathbf{X}},\mathbf{X})}{1+{\delta }_{\mathrm{sin}}}$, with ${\mathbf{K}}_{\mathrm{cos}}(\hat{\mathbf{X}},\mathbf{X}),{\mathbf{K}}_{\mathrm{sin}}(\hat{\mathbf{X}},\mathbf{X})\in {\mathbb{R}}^{\hat{n}\times n}$ and ${\mathbf{K}}_{\mathrm{cos}}(\hat{\mathbf{X}},\hat{\mathbf{X}}),{\mathbf{K}}_{\mathrm{sin}}(\hat{\mathbf{X}},\hat{\mathbf{X}})\in {\mathbb{R}}^{\hat{n}\times \hat{n}}$ defined as in (8).

Proof. See appendix C. □

Similar to theorem 2 on ${\bar{E}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}$, here the expression for ${\bar{E}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}$ is also given as the sum of first- and second-order corrections. To see this, one can confirm, by taking $(\hat{\mathbf{X}},\hat{\mathbf{y}})=(\mathbf{X},\mathbf{y})$, that the first term in ${\bar{E}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}$ becomes

$$\begin{equation*}\frac{1}{\hat{n}}{\Vert}\hat{\mathbf{y}}-\frac{N}{n}\hat{\mathbf{\Phi }}\bar{\mathbf{Q}}\mathbf{y}{{\Vert}}^{2}=\frac{1}{n}{\Vert}\mathbf{y}-\frac{N}{n}\mathbf{\Phi }\bar{\mathbf{Q}}\mathbf{y}{{\Vert}}^{2}=\frac{{\lambda }^{2}}{n}{\Vert}\bar{\mathbf{Q}}\mathbf{y}{{\Vert}}^{2}\end{equation*}$$

and is equal to the first term in ${\bar{E}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}$, where we used the fact that $\frac{N}{n}\mathbf{\Phi }\bar{\mathbf{Q}}={\mathbf{I}}_{n}-\lambda \bar{\mathbf{Q}}$. The same also holds for the second term, so that one obtains ${\bar{E}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}={\bar{E}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}$, with $(\hat{\mathbf{X}},\hat{\mathbf{y}})=(\mathbf{X},\mathbf{y})$, as expected. From this perspective, theorem 3 can be seen as an extension of theorem 2, with the 'interaction' between training and test data (e.g. test-versus-test ${\mathbf{K}}_{\sigma }(\hat{\mathbf{X}},\hat{\mathbf{X}})$ and test-versus-train ${\mathbf{K}}_{\sigma }(\hat{\mathbf{X}},\mathbf{X})$ interaction matrices) summarized in the scalar parameter Θ_σ defined in (14), for σ ∈ { cos,  sin }.

By taking N/n → ∞, we have that $\bar{\mathbf{Q}}\sim \frac{n}{N}{\mathbf{K}}^{-1}$, Θ_σ ∼ N⁻¹, Ω → I₂, and consequently

$$\begin{equation*}\underset{N/n\to \infty }{\mathrm{lim}}\enspace {\bar{E}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}=\frac{1}{\hat{n}}{\Vert}\hat{\mathbf{y}}-\mathbf{K}(\hat{\mathbf{X}},\mathbf{X}){\mathbf{K}}_{\mathbf{X}}^{-1}\mathbf{y}{{\Vert}}^{2}.\end{equation*}$$

This is the test MSE of classical Gaussian kernel regression, with $\mathbf{K}(\hat{\mathbf{X}},\mathbf{X})\equiv {\mathbf{K}}_{\mathrm{cos}}(\hat{\mathbf{X}},\mathbf{X})+{\mathbf{K}}_{\mathrm{sin}}(\hat{\mathbf{X}},\mathbf{X})\in {\mathbb{R}}^{\hat{n}\times n}$ the test-versus-train Gaussian kernel matrix. As opposed to the training MSE discussed in remark 4, here ${\bar{E}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}$ generally has a non-zero limit (that is, however, independent of λ) as N/n → ∞.

The RFF Gram matrix ${\mathbf{\Sigma }}_{\mathbf{X}}^{\mathsf{\text{T}}}{\mathbf{\Sigma }}_{\mathbf{X}}/N$ is not close to the classical Gaussian kernel matrix K_X in the large n, p, N regime; and, as a consequence, its resolvent Q, as well the training and test MSE, E_train and E_test (that are functions of Q), behave quite differently from the Gaussian kernel predictions. As already discussed in remark 2 after theorem 1, for λ > 0, the pair (δ_cos, δ_sin) characterizes the correction when considering n, p, N all large, compared to the large-N-only asymptotic behavior:

$$\begin{equation}{\delta }_{\mathrm{cos}}=\frac{1}{n}\enspace \mathrm{t}\mathrm{r}\enspace {\mathbf{K}}_{\mathrm{cos}}\bar{\mathbf{Q}},{\delta }_{\mathrm{sin}}=\frac{1}{n}\enspace \mathrm{t}\mathrm{r}\enspace {\mathbf{K}}_{\mathrm{sin}}\bar{\mathbf{Q}},\quad \bar{\mathbf{Q}}={\left(\frac{N}{n}\left(\frac{{\mathbf{K}}_{\mathrm{cos}}}{1+{\delta }_{\mathrm{cos}}}+\frac{{\mathbf{K}}_{\mathrm{sin}}}{1+{\delta }_{\mathrm{sin}}}\right)+\lambda {\mathbf{I}}_{n}\right)}^{-1}.\end{equation} \tag{ 15 }$$

To start, figure 2 compares the training MSEs of RFF ridge regression to the predictions from Gaussian kernel regression and to the predictions from our theorem 2, on the popular MNIST data set [47]. Observe that there is a huge gap between empirical training errors and the Gaussian kernel predictions, especially when N/n < 1, while our theory consistently fits empirical observations almost perfectly.

Next, from (15) we know that both δ_cos and δ_sin are decreasing functions of λ. (See lemma 7 in appendix D for a proof of this fact.) Figure 3 shows that: (i) over a range of different N/n, both δ_cos and δ_sin decrease monotonically as λ increases; (ii) the behavior for N/n < 1, which is decreasing from an initial value of δ ≫ 1, is very different from the behavior for N/n ≳ 1, where an initially flat region is observed for small values of λ and we have δ < 1 for all values of λ; and (iii) the impact of regularization λ becomes less significant as the ratio N/n becomes large. This is in accordance with the limiting behavior of $\bar{\mathbf{Q}}\simeq \frac{n}{N}{\mathbf{K}}_{\mathbf{X}}^{-1}$ in remark 2 that is independent of λ as N/n → ∞.

Note also that, while δ_cos and δ_sin can be geometrically interpreted as a sort of weighted 'angle' between different kernel matrices, and therefore one might expect to have δ ∈ [0, 1], this is not the case for the leftmost plot with N/n = 1/4. There, for small values of λ (say λ ≲ 0.1), both δ_cos and δ_sin scale like λ⁻¹, while they are observed to saturate to a fixed O(1) value for N/n = 1, 4, 16. This corresponds to two different phases of learning in the 'ridgeless' λ → 0 case. As we shall see in more detail later in section 4.1, depending on whether we are in the 'under-parameterized' (2N < n) or the 'over-parameterized' (2N > n) regime, the system behaves fundamentally differently.

Both δ_cos and δ_sin in (15) are decreasing functions of N, as depicted in figure 4. (See lemma 6 in appendix D for a proof.) More importantly, figure 4 also illustrates that δ_cos and δ_sin exhibit qualitatively different behavior: for λ not too small (λ = 1 or 10), we observe a rather 'smooth' behavior, as a function of the ratio N/n, and they both decrease smoothly, as N/n grows large. However, for λ relatively small (λ = 10⁻³ and 10⁻⁷), we observe a sharp 'phase transition' on two sides of the interpolation threshold 2N = n. (Note that the scale of the y-axis is very different in different subfigures.) More precisely, in the leftmost plot with λ = 10⁻⁷, the values of δ_cos and δ_sin 'jumps' from order O(1) (when 2N > n) to much higher values of the order of λ⁻¹ (when 2N < n). A similar behavior is also observed for λ = 10⁻³.

As a consequence of this phase transition, different behaviors are expected for training and test MSEs in the 2N < n and 2N > n regime. Figure 5 depicts the empirical and theoretical test MSEs with different regularization penalty λ. In particular, for λ = 10⁻⁷ and λ = 10⁻³, a double descent behavior is observed, with a singularity at 2N = n, while for larger values of λ (λ = 0.2, 10), a smoother and monotonically decreasing curve for test error is observed, as a function of N/n. Figure 5 also illustrates that: (i) for a fixed regularization λ > 0, the minimum test error is always obtained in the over-parameterization 2N > n regime; and (ii) the global optimal design (over N and λ) is achieved by highly over-parametrized system with a (problem-dependent) non-vanishing λ. This is in accordance with the observations in [12] for Gaussian data.

Remark 5. (On ridge regularization). Performing ridge regularization (with λ as a control parameter) is known to help alleviate the sharp performance drop around 2N = n [12, 18]. Our theorem 3 can serve as a convenient alternative to evaluate the effect of small λ around 2N = n, as well as to determine an optimal λ, for not-too-small n, p, N. In the setup of figure 5, a grid search can be used to find the regularization that minimizes ${\bar{E}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}$. For this choice of λ (λ_opt ≈ 0.2), no singular peak at 2N = n is observed.

Remark 6. (Double descent as a consequence of phase transition). While the double descent phenomenon has received considerable attention recently, our analysis makes it clear that in this model (and presumably many others) it is a natural consequence of the phase transition between two qualitatively different phases of learning [7].

We chose to present our theoretical results in section 2 (theorems 1–3) in the same form, regardless of whether 2N > n or 2N < n. This comes at the cost of requiring a strictly positive ridge regularization λ > 0, as n, p, N → ∞. As discussed in section 3, for small values of λ, depending on the sign of 2N − n, we observe totally different behaviors for (δ_cos, δ_sin) and thus for the key resolvent $\bar{\mathbf{Q}}$. As a matter of fact, for λ = 0 and 2N < n, the (random) resolvent Q(λ = 0) in (7) is simply undefined, as it involves inverting a singular matrix ${\mathbf{\Sigma }}_{\mathbf{X}}^{\mathsf{\text{T}}}{\mathbf{\Sigma }}_{\mathbf{X}}\in {\mathbb{R}}^{n\times n}$ that is of rank at most 2N < n. As a consequence, we expect to see $\bar{\mathbf{Q}}\sim {\lambda }^{-1}$ as λ → 0 for 2N < n, while for 2N > n this is not the case.

These two phases of learning can be theoretically justified by considering the ridgeless λ → 0 limit in theorem 1, with the unified variables γ_cos and γ_sin introduced below.

• (a)  
  For 2N < n and λ → 0, we obtain

  $$\begin{equation}\begin{cases}\lambda {\delta }_{\mathrm{cos}}\to {\gamma }_{\mathrm{cos}}\equiv \frac{1}{n}\enspace \mathrm{t}\mathrm{r}\enspace {\mathbf{K}}_{\mathrm{cos}}{\left(\frac{N}{n}\left(\frac{{\mathbf{K}}_{\mathrm{cos}}}{{\gamma }_{\mathrm{cos}}}+\frac{{\mathbf{K}}_{\mathrm{sin}}}{{\gamma }_{\mathrm{sin}}}\right)+{\mathbf{I}}_{n}\right)}^{-1}\quad \hfill \\ \lambda {\delta }_{\mathrm{sin}}\to {\gamma }_{\mathrm{sin}}\equiv \frac{1}{n}\enspace \mathrm{t}\mathrm{r}\enspace {\mathbf{K}}_{\mathrm{sin}}{\left(\frac{N}{n}\left(\frac{{\mathbf{K}}_{\mathrm{cos}}}{{\gamma }_{\mathrm{cos}}}+\frac{{\mathbf{K}}_{\mathrm{sin}}}{{\gamma }_{\mathrm{sin}}}\right)+{\mathbf{I}}_{n}\right)}^{-1}\quad \hfill \end{cases},\end{equation} \tag{ 16 }$$

  in such as way that δ_cos, δ_sin and $\bar{\mathbf{Q}}$ scale like λ⁻¹. We have in particular $\mathbb{E}[\lambda \mathbf{Q}]\sim \lambda \bar{\mathbf{Q}}\sim {\left(\frac{N}{n}\left(\frac{{\mathbf{K}}_{\mathrm{cos}}}{{\gamma }_{\mathrm{cos}}}+\frac{{\mathbf{K}}_{\mathrm{sin}}}{{\gamma }_{\mathrm{sin}}}\right)+{\mathbf{I}}_{n}\right)}^{-1}$ with (γ_cos, γ_sin) of order O(1).
• (b)  
  For 2N > n and λ → 0, we obtain

  $$\begin{equation}\begin{cases}{\delta }_{\mathrm{cos}}\to {\gamma }_{\mathrm{cos}}=\frac{1}{N}\enspace \mathrm{t}\mathrm{r}\enspace {\mathbf{K}}_{\mathrm{cos}}{\left(\frac{{\mathbf{K}}_{\mathrm{cos}}}{1+{\gamma }_{\mathrm{cos}}}+\frac{{\mathbf{K}}_{\mathrm{sin}}}{1+{\gamma }_{\mathrm{sin}}}\right)}^{-1}\quad \hfill \\ {\delta }_{\mathrm{sin}}\to {\gamma }_{\mathrm{cos}}=\frac{1}{N}\enspace \mathrm{t}\mathrm{r}\enspace {\mathbf{K}}_{\mathrm{sin}}{\left(\frac{{\mathbf{K}}_{\mathrm{cos}}}{1+{\gamma }_{\mathrm{cos}}}+\frac{{\mathbf{K}}_{\mathrm{sin}}}{1+{\gamma }_{\mathrm{sin}}}\right)}^{-1}\quad \hfill \end{cases},\end{equation} \tag{ 17 }$$

  by taking directly λ → 0 in theorem 1.

As a consequence, in the ridgeless limit λ → 0, theorem 1 exhibits the following two learning phases:

• (a)  
  Under-parameterized phase: with 2N < n. Here, Q is not well defined (indeed Q ∼ λ⁻¹) and one must consider instead the properly scaled γ_cos, γ_sin and $\lambda \bar{\mathbf{Q}}$ in (16). Like δ_cos and δ_sin, γ_cos and γ_sin also decrease as N/n grows large. In particular, one has ${\gamma }_{\mathrm{cos}},{\gamma }_{\mathrm{sin}},{\Vert}\lambda \bar{\mathbf{Q}}{\Vert}\to 0$ as 2N − n ↑ 0.
• (b)  
  Over-parameterized phase: with 2N > n. Here, one can consider δ_cos, δ_sin and ${\Vert}\bar{\mathbf{Q}}{\Vert}$. One has particularly that ${\delta }_{\mathrm{cos}},{\delta }_{\mathrm{sin}},{\Vert}\bar{\mathbf{Q}}{\Vert}\to \infty$ as 2N − n ↓ 0 and tend to zero as N/n → ∞.

With this discussion on the two phases of learning, we now understand why:

• in the leftmost plot of figure 3 with 2N < n, δ_cos and δ_sin behave rather differently from other plots and approximately scale as λ⁻¹ for small values of λ; and
• in the first and second leftmost plots of figure 4, a 'jump' in the values of δ occurs at the transition point 2N = n, and the δ's are numerically of the same order of λ⁻¹ for 2N < n.

To characterize the phase transition from (16) and (17) in the λ → 0 setting, we consider the scaled variables

$$\begin{equation}\begin{cases}{\gamma }_{\sigma }=\lambda {\delta }_{\sigma }\quad \hfill & \text{for}\enspace 2\enspace N< n\hfill \\ {\gamma }_{\sigma }={\delta }_{\sigma }\quad \hfill & \text{for}\enspace 2\enspace N > n\hfill \end{cases},\quad \sigma \in \left\{\mathrm{cos},\mathrm{sin}\right\}.\end{equation} \tag{ 18 }$$

An advantage of using these scaled variables is that they are of order O(1) as n, p, N → ∞ and λ → 0. The behavior of (γ_cos, γ_sin) is reported in figure 6, in the same setting as figure 4. Observe the sharp transition between the 2N < n and 2N > n regime, in particular for λ = 10⁻⁷ and λ = 10⁻³, and that this transition is smoothed out for λ = 1. (A 'transition' is also seen for λ = 10, but this is potentially misleading. It is true that γ_cos and γ_sin do change in this way, as a function of N/n, but unless λ ≈ 0, these quantities are not solutions of the aforementioned fixed point equations.)

Zoom In Zoom Out Reset image size

On account of these two different phases of learning (under- and over-parameterized, in (16) and (17), respectively) and the sharp transition of (γ_cos, γ_sin) in figure 6, it is not surprising to observe a 'singular' behavior at 2N = n, when no regularization is applied. We next examine the asymptotic training and test errors in more detail.

Asymptotic training MSE as λ → 0. In the under-parameterized regime with 2N < n, combining (16) we have that both $\lambda \bar{\mathbf{Q}}$ and $\frac{\bar{\mathbf{Q}}}{1+{\delta }_{\sigma }}\sim \frac{\lambda \bar{\mathbf{Q}}}{{\gamma }_{\sigma }},\sigma \in \left\{\mathrm{cos},\mathrm{sin}\right\}$ are well-behaved and are generally not zero. As a consequence, by theorem 2, the asymptotic training error ${\bar{E}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}$ tends to a nonzero limit as λ → 0, measuring the residual information in the training set that is not captured by the regressor $\boldsymbol{\beta }\in {\mathbb{R}}^{2N}$. As 2N − n ↑ 0, we have γ_cos, γ_sin → 0 and ${\Vert}\lambda \bar{\mathbf{Q}}{\Vert}\to 0$ so that ${\bar{E}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}\to 0$ and β interpolates the entire training set. On the other hand, in the over-parameterized 2N > n regime, one always has ${\bar{E}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}=0$. This particularly implies the training error is 'continuous' around the point 2N = n.

Asymptotic test MSE as λ → 0. Again, in the under-parameterized regime with 2N < n, now consider the more involved asymptotic test error in theorem 3. In particular, we will focus here on the case $\hat{\mathbf{X}}\ne \mathbf{X}$ (or, more precisely, they are sufficiently different from each other in such a way that ${\Vert}\mathbf{X}-\hat{\mathbf{X}}{\Vert}\to \hspace{-12pt}/\;\;\;0$ as n, p, N → ∞ and λ → 0; see further discussion below in section 4.2) so that ${\mathbf{K}}_{\sigma }(\mathbf{X},\mathbf{X})\ne {\mathbf{K}}_{\sigma }(\hat{\mathbf{X}},\mathbf{X})$ and $\frac{N}{n}\hat{\mathbf{\Phi }}\bar{\mathbf{Q}}\ne {\mathbf{I}}_{n}-\lambda \bar{\mathbf{Q}}$. In this case, the two-by-two matrix Ω in ${\bar{E}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}$ diverges to infinity at 2N = n in the λ → 0 limit. (Indeed, the determinant det(Ω⁻¹) scales as λ, per lemma 5.) As a consequence, we have ${\bar{E}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}\to \infty$ as 2N → n, resulting in a sharp deterioration of the test performance around 2N = n. (Of course, this holds if no additional regularization is applied as discussed in remark 5.) It is also interesting to note that, while Ω also appears in ${\bar{E}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}$, we still obtain (asymptotically) zero training MSE at 2N = n, despite the divergence of Ω, again due to the prefactor λ² in ${\bar{E}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}$. If λ ≳ 1, then det(Ω⁻¹) exhibits much more regular properties (figure 7), as one would expect.

Zoom In Zoom Out Reset image size

Continuing our discussion of the RFF performance in the large n, p, N limit, we can see that the (asymptotic) test error behaves entirely differently, depending on whether $\hat{\mathbf{X}}$ is 'close to' X or not. For $\hat{\mathbf{X}}=\mathbf{X}$, one has ${\bar{E}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}={\bar{E}}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}$ that decreases monotonically as N grows large; while for $\hat{\mathbf{X}}$ 'sufficiently' different from X, ${\bar{E}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}$ diverges to infinity at 2N = n. To have a more quantitative assessment of the influence of training-test data similarity on the test error, we consider the special case $\hat{n}=n$ and $\hat{\mathbf{y}}=\mathbf{y}$. In this case, it follows from theorem 3 that

$$\begin{align*}\hfill {{\Theta}}_{\sigma }& =\frac{1}{N}\enspace \mathrm{t}\mathrm{r}({\mathbf{K}}_{\sigma }+{\mathbf{K}}_{\sigma }(\hat{\mathbf{X}},\hat{\mathbf{X}})-2{\mathbf{K}}_{\sigma }(\hat{\mathbf{X}},\mathbf{X}))+\frac{2}{n}\enspace \mathrm{t}\mathrm{r}\enspace \bar{\mathbf{Q}}{\Delta}{\mathbf{\Phi }}^{\mathsf{\text{T}}}{\Delta}{\mathbf{K}}_{\sigma }\hfill \\ \hfill & \quad +\frac{N}{n}\frac{1}{n}\enspace \mathrm{t}\mathrm{r}\enspace \bar{\mathbf{Q}}{\Delta}{\mathbf{\Phi }}^{\mathsf{\text{T}}}{\Delta}\mathbf{\Phi }\bar{\mathbf{Q}}{\mathbf{K}}_{\sigma }+\frac{n}{N}\frac{{\lambda }^{2}}{n}\enspace \mathrm{t}\mathrm{r}\enspace \bar{\mathbf{Q}}{\mathbf{K}}_{\sigma }\bar{\mathbf{Q}}-\frac{2\lambda }{N}\enspace \mathrm{t}\mathrm{r}\enspace \bar{\mathbf{Q}}{\Delta}{\mathbf{K}}_{\sigma }-\frac{2\lambda }{n}\enspace \mathrm{t}\mathrm{r}\enspace \bar{\mathbf{Q}}{\Delta}{\mathbf{\Phi }}^{\mathsf{\text{T}}}\bar{\mathbf{Q}}{\mathbf{K}}_{\sigma },\hfill \end{align*}$$

for σ ∈ { cos,  sin }, ${\Delta}{\mathbf{K}}_{\sigma }={\mathbf{K}}_{\sigma }-{\mathbf{K}}_{\sigma }(\hat{\mathbf{X}},\mathbf{X})$ and ${\Delta}\mathbf{\Phi }\equiv \hat{\mathbf{\Phi }}-\mathbf{\Phi }$. Since in the ridgeless λ → 0 limit the matrix Ω scale as λ⁻¹ (see figure 7), one must have Θ_σ scaling as λ so that ${\bar{E}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}$ does not diverge at 2N = n as λ → 0. One example is the case where the test data is a small (additive) perturbation of the training data such that, in the kernel feature space

$$\begin{equation*}{\mathbf{K}}_{\sigma }-{\mathbf{K}}_{\sigma }(\hat{\mathbf{X}},\mathbf{X})=\lambda {\mathbf{\Xi }}_{\sigma },\quad {\mathbf{K}}_{\sigma }(\hat{\mathbf{X}},\hat{\mathbf{X}})-{\mathbf{K}}_{\sigma }(\hat{\mathbf{X}},\mathbf{X})=\lambda {\hat{\mathbf{\Xi }}}_{\sigma }\end{equation*}$$

for ${\mathbf{\Xi }}_{\sigma },{\hat{\mathbf{\Xi }}}_{\sigma }\in {\mathbb{R}}^{n\times n}$ of bounded spectral norms. In this setting, we have ${{\Theta}}_{\sigma }=\frac{\lambda }{N}\enspace \mathrm{t}\mathrm{r}({\mathbf{\Xi }}_{\sigma }+{\hat{\mathbf{\Xi }}}_{\sigma })+O({\lambda }^{2})$ so that the asymptotic test error does not diverge to infinity at 2N = n as λ → 0. This is supported by figure 8, where the test data are generated by adding Gaussian white noise of variance σ² to the training data, i.e. ${\hat{\mathbf{x}}}_{i}={\mathbf{x}}_{i}+\sigma {\boldsymbol{\varepsilon }}_{i}$, for independent ${\boldsymbol{\varepsilon }}_{i}\sim \mathcal{N}(\mathbf{0},{\mathbf{I}}_{p}/p)$. In figure 8, we observe that (i) below the threshold σ² = λ, test error coincides with the training error and both are close to zero; and (ii) as soon as σ² ≳ λ, the test error diverges from the training error and grows large (but linearly in σ²) as the noise level increases. Note also from the two rightmost plots of figure 8 that, the training-to-test 'transition' at σ² ∼ λ is sharp only for relatively small values of λ, as predicted by our theory.

Zoom In Zoom Out Reset image size

So far, we have presented results in detail for one particular real-world data set, but we have extensive empirical results demonstrating that similar conclusions hold more broadly. As an example of these additional results, here we present a numerical evaluation of our results on several other real-world image data sets. We consider the classification task on another two MNIST-like data sets composed of 28 × 28 grayscale images: the Fashion-MNIST [48] and the Kannada-MNIST [49] data sets. Each image is represented as a p = 784-dimensional vector and the output targets $\mathbf{y},\hat{\mathbf{y}}$ are taken to have −1, +1 entries depending on the image class. As a consequence, both the training and test MSEs in (5) are approximately 1 for N = 0 and significantly small λ, as observed in figures 5 and 11 below. For each data set, images were jointly centered and scaled so to fall close to the setting of assumption 1 on X and $\hat{\mathbf{X}}$.

In figure 9, we compare the empirical training and test errors with their limiting behaviors derived in theorems 2 and 3, as a function of the penalty parameter λ, on a training set of size n = 1024 (512 images from class 5 and 512 images from class 6) with feature dimension N = 256 and N = 512, on both data sets. A close fit between theory and practice is observed, for moderately large values of n, p, N, demonstrating a wide practical applicability of the proposed asymptotic analyses, particularly compared to the (limiting) Gaussian kernel predictions per figure 2.

Zoom In Zoom Out Reset image size

In figure 10, we report the behavior of the pair (δ_cos, δ_sin) for small values of λ = 10⁻⁷ and 10⁻³. Similar to the two leftmost plots in figure 4 for MNIST, a jump from the under-to over-parameterized regime occurs at the interpolation threshold 2N = n, in both Fashion- and Kannada-MNIST data sets, clearly indicating the two phases of learning and the phase transition between them.

Zoom In Zoom Out Reset image size

In figure 11, we report the empirical and theoretical test errors as a function of the ratio N/n, on a training test of size n = 500 (250 images from class 8 and 250 images from class 9), by varying the feature dimension N. An exceedingly small regularization λ = 10⁻⁷ is applied to mimic the 'ridgeless' limiting behavior as λ → 0. On both data sets, the corresponding double descent curve is observed where the test errors go down and up, with a singular peak around 2N = n, and then go down monotonically as N continues to increase when 2N > n.

Zoom In Zoom Out Reset image size