Harmless interpolation of noisy data in regression

TL;DR

This paper analyzes how overparameterized linear regression models can interpolate noisy data without harming generalization, showing that increasing features can reduce mean-squared error and proposing a hybrid scheme for sparse models.

Contribution

It characterizes the generalization error of interpolating solutions in noisy linear regression and introduces a hybrid method to improve performance in sparse models.

Findings

Overparameterization can lead to harmless noise interpolation.

Generalization error decreases as the number of features increases.

A hybrid interpolating scheme achieves near-optimal MSE in sparse models.

Abstract

A continuing mystery in understanding the empirical success of deep neural networks is their ability to achieve zero training error and generalize well, even when the training data is noisy and there are more parameters than data points. We investigate this overparameterized regime in linear regression, where all solutions that minimize training error interpolate the data, including noise. We characterize the fundamental generalization (mean-squared) error of any interpolating solution in the presence of noise, and show that this error decays to zero with the number of features. Thus, overparameterization can be explicitly beneficial in ensuring harmless interpolation of noise. We discuss two root causes for poor generalization that are complementary in nature -- signal "bleeding" into a large number of alias features, and overfitting of noise by parsimonious feature selectors. For the sparse linear model with noise, we provide a hybrid interpolating scheme that mitigates both these issues and achieves order-optimal MSE over all possible interpolating solutions.

Figures

Tables

Table 1. Table 1: Corollaries of Proposition  1 for a variety of initial estimators for the ( k , σ ) 𝑘 𝜎 (k,\sigma) -sparse linear model.

Estimator ${\widehat{𝜶}}_1$	Side information	Test MSE of hybrid interpolator
Optimal (e.g. SLOPE)	${\sigma }^2$	$𝒪\left({\sigma }^2\cdot \frac{k\ln \left(\frac{d}{k}\right)}{n}+{\sigma }^2\cdot \frac{n}{d}\right)+{\sigma }^2$
Lagrangian Lasso	${\sigma }^2$	$𝒪\left({\sigma }^2\cdot \frac{k\ln d}{n}+{\sigma }^2\cdot \frac{n}{d}\right)+{\sigma }^2$
OMP upto $k$ steps	$k$	$𝒪\left({\sigma }^2\cdot \frac{k\ln d}{n}+{\sigma }^2\cdot \frac{n}{d}\right)+{\sigma }^2$
OMP upto $k_0>k$ steps	Upper bound on $k$	$𝒪\left({\sigma }^2\cdot \frac{k_0\ln d}{n}+{\sigma }^2\cdot \frac{n}{d}\right)+{\sigma }^2$
Square-root-Lasso	None	$𝒪\left({\sigma }^2\cdot \frac{k\ln d}{n}+{\sigma }^2\cdot \frac{n}{d}\right)+{\sigma }^2$

Equations

$${\mathbf{A}}_{\mathsf{train}}:={\left[\begin{array}{ccc}\mathbf{a}({\mathbf{X}}_1{)}^{\top }& \mathbf{a}({\mathbf{X}}_2{)}^{\top }& \dots \mathbf{a}({\mathbf{X}}_n{)}^{\top }\end{array}\right]}^{\top }$$

$${\mathbf{A}}_{\mathsf{train}}\mathbf{\alpha }={\mathbf{Y}}_{\mathsf{train}}$$

$${\mathcal{E}}_{\mathsf{test}}(\hat{\mathbf{\alpha }}):=\mathbb{E}[(Y-⟨\mathbf{a}(\mathbf{X}),\hat{\mathbf{\alpha }}⟩{)}^2]-{\sigma }^2$$

$${\hat{\mathbf{\alpha }}}_{\mathsf{ideal}}:=\arg \min {\mathcal{E}}_{\mathsf{test}}(\mathbf{\alpha })$$

$$\text{ subject to }{\mathbf{A}}_{\mathsf{train}}\mathbf{\alpha }={\mathbf{Y}}_{\mathsf{train}}.$$

$${\mathcal{E}}_{\mathsf{test}}^{\ast }={\mathbf{W}}_{\mathsf{train}}^{\top }({\mathbf{B}}_{\mathsf{train}}{\mathbf{B}}_{\mathsf{train}}^{\top }{)}^{-1}{\mathbf{W}}_{\mathsf{train}}.$$

$$\mathbf{a}(X)=\left[\begin{array}{ccccc}1& e^{2\pi iX}& e^{2\pi (2i)X}& \dots & e^{2\pi ((d-1)i))X}\end{array}\right]\in {\mathbb{C}}^d.$$

$$\mathbf{a}(X)=\left[\begin{array}{cccc}1& X& X^2& \dots X^{d-1}.\end{array}\right]$$

$$\mathbf{a}(X)=\left[\begin{array}{ccc}{p}_0(X)& p_1(X)& \dots p_{d-1}(X),\end{array}\right]$$

$${\mathcal{E}}_{\mathsf{test}}^{\ast }\ge \left(\frac{n(1-\delta )}{(C\sqrt{d\ln n}+\sqrt{n}{)}^2}+1\right){\sigma }^2=\omega \left(\frac{{\sigma }^{2}n}{d\ln n}\right)$$

$${\mathcal{E}}_{\mathsf{test}}^{\ast }\ge \left(\frac{n(1-\delta )}{(C_K\sqrt{d}+\sqrt{n}{)}^2}+1\right){\sigma }^2=\omega \left(\frac{{\sigma }^{2}n}{d}\right)$$

$${\mathcal{E}}_{\mathsf{test}}^{\ast }\ge \left(\frac{n(1-\delta )}{(\sqrt{d}+2\sqrt{n}{)}^2}+1\right){\sigma }^2=\omega \left(\frac{{\sigma }^{2}n}{d}\right)$$

$${\mathcal{E}}_{\mathsf{test}}^{\ast }\le \left(\frac{n(1+\delta )}{(\sqrt{d}-2\sqrt{n}{)}^2}+1\right){\sigma }^2=\mathcal{O}\left(\frac{{\sigma }^{2}n}{d}\right)$$

$${\mathcal{E}}_{\mathsf{test}}^{\ast }\le \left(\frac{4C_K^{2}n(1+\delta )}{(\sqrt{d}-\sqrt{n-1}{)}^2}+1\right){\sigma }^2=\mathcal{O}\left(\frac{{\sigma }^{2}n}{d}\right)$$

$${\mathcal{E}}_{\mathsf{test}}(\hat{\mathbf{\alpha }})$$

$$=\mathbb{E}[(⟨\mathbf{a}(\mathbf{X}),\hat{\mathbf{\alpha }}-{\mathbf{\alpha }}^{\ast }⟩{)}^2]+\mathbb{E}[W^2]$$

$$=\mathbb{E}[(\hat{\mathbf{\alpha }}-{\mathbf{\alpha }}^{\ast }{)}^{\top }\mathbf{a}(\mathbf{X})\mathbf{a}(\mathbf{X}{)}^{\top }(\hat{\mathbf{\alpha }}-{\mathbf{\alpha }}^{\ast })]+{\sigma }^2-{\sigma }^2$$

$$=(\hat{\mathbf{\alpha }}-{\mathbf{\alpha }}^{\ast }{)}^{\top }\mathbf{\Sigma }(\hat{\mathbf{\alpha }}-{\mathbf{\alpha }}^{\ast })+{\sigma }^2-{\sigma }^2$$

$$=\Vert {\mathbf{\Sigma }}^{1/2}(\hat{\mathbf{\alpha }}-{\mathbf{\alpha }}^{\ast }){\Vert }_2^2+{\sigma }^2-{\sigma }^2$$

$$=\Vert {\mathbf{\Sigma }}^{1/2}(\hat{\mathbf{\alpha }}-{\mathbf{\alpha }}^{\ast }){\Vert }_2^2.$$

$${\hat{\mathbf{\alpha }}}_{\mathsf{ideal}}$$

$${\mathbf{A}}_{\mathsf{train}}(\mathbf{\alpha }-{\mathbf{\alpha }}^{\ast })$$

$$⟹{\mathbf{A}}_{\mathsf{train}}{\mathbf{\Sigma }}^{-1/2}{\mathbf{\Sigma }}^{1/2}(\mathbf{\alpha }-{\mathbf{\alpha }}^{\ast })$$

$${\mathbf{\Sigma }}^{1/2}({\hat{\mathbf{\alpha }}}_{\mathsf{ideal}}-{\mathbf{\alpha }}^{\ast })={\mathbf{B}}_{\mathsf{train}}^{\top }({\mathbf{B}}_{\mathsf{train}}{\mathbf{B}}_{\mathsf{train}}^{\top }{)}^{-1}{\mathbf{W}}_{\mathsf{train}}$$

$${\mathcal{E}}_{\mathsf{test}}^{\ast }$$

$$={\mathbf{W}}_{\mathsf{train}}^{\top }({\mathbf{B}}_{\mathsf{train}}{\mathbf{B}}_{\mathsf{train}}^{\top }{)}^{-1}{\mathbf{W}}_{\mathsf{train}}$$

$$\sqrt{d}-\sqrt{n}-t\le {\sigma }_{min}(\mathbf{B})\le {\sigma }_{max}(\mathbf{B})\le \sqrt{d}+\sqrt{n}+t$$

$${\mathcal{E}}_{\mathsf{test}}^{\ast }$$

$$\Vert {\mathbf{W}}_{\mathsf{train}}{\Vert }_2^2\ge n{\sigma }^2(1-\delta )$$

$${\mathcal{E}}_{\mathsf{test}}^{\ast }$$

$$\Vert {\mathbf{W}}_{\mathsf{train}}{\Vert }_2^2\le n{\sigma }^2(1+\delta )$$

$$Y=1+W$$

Full text

1 Introduction

In statistical machine learning, we have high-dimensional data in the form of $n$ covariate-response pairs $(\mathbf{X}_{i}\in\mathbb{R}^{d},Y_{i})_{i=1}^{n}$. When training parametric models for functions fitting covariate to response, the traditional wisdom [1] is to select function classes with a number of parameters $d<n$. In classification (when the labels are discrete), the scaling of the test error with respect to $n$ is determined by the VC-dimension [2]/Rademacher complexity [3] of the function class, which in the worst case increases with the number of parameters $d$. In regression (when the labels are continuous), the mean-squared error (MSE) of the least-squares estimator scales as the condition number of the data matrix, which is reasonable for smaller ratios of $d/n$ but increases astronomically as $d$ approaches $n$. The qualitative fear is the same: if the function class is too complex, it starts to overfit noise and generalizes poorly to unseen test data.

This wisdom has been challenged by the recent advent of deeper and deeper neural networks. In particular, a thought-provoking paper [4] noted that several deep neural networks generalize well despite achieving zero or close to zero training error, and being so expressive that they even have the ability to fit pure noise. As they put it, “understanding deep learning requires rethinking generalization". How can we reconcile the fact that good interpolating solutions exist with the classical bias-variance tradeoff?

These phenomena are being actively investigated for neural networks and kernel methods. A natural simpler setting to investigate is overparameterized linear regression on noisy data. The overparameterized (high-dimensional) linear model has the advantage of being classically studied under estimation procedures that incorporate some kind of regularization111The most commonly known are Tikhonov (ridge) regularization and $\ell_{1}$-regularization, but we refer the interested reader to Miller’s book [5] for a comprehensive list.. Such procedures typically avoid fitting noise, thus incurring some non-zero training error. We are now interested in whether we can achieve any success with solutions that interpolate the training data. In particular, can we ever interpolate the noise and retain favorable generalization guarantees? If so, which interpolators should we use, and which should we avoid?

In this paper, we provide constructive answers to the above questions for overparameterized linear regression using elementary machinery. Our contributions are as follows:

We give a fundamental limit (Theorem 1, Corollaries 1 and 2) for the excess MSE of any interpolating solution in the presence of noise, and show that it converges to [math] as the number of features goes to infinity for feature families satisfying mild conditions. 2. 2.

We provide a Fourier-theoretic interpretation of concurrent analyses [6, 7, 8, 9, 10] of the minimum $\ell_{2}$-norm interpolator. 3. 3.

We show (Theorem 2) that parsimonious interpolators (like the $\ell_{1}$-minimizing interpolator and its relatives) suffer the complementary problem of overfitting pure noise. 4. 4.

We construct two-step hybrid interpolators that successfully recover signal and harmlessly fit noise, achieving the order-optimal rate of test MSE among all interpolators (Proposition 1 and all its corollaries).

1.1 Related work

We discuss prior work in three categories: a) overparameterization in deep neural networks, b) interpolation of high-dimensional data using kernels, and c) high-dimensional linear regression. We then recap work on overparameterized linear regression that is concurrent to ours.

1.1.1 Recent interest in overparameterization

Conventional statistical wisdom is that using more parameters in one’s model than data points leads to poor generalization. This wisdom is corroborated in theory by worst-case generalization bounds on such overparameterized models following from VC-theory in classification [2] and ill-conditioning in least-squares regression [5]. It is, however, contradicted in practice by the notable recent trend of empirically successful overparameterized deep neural networks. For example, the commonly used CIFAR-$10$ dataset contains $60000$ images, but the number of parameters in all the neural networks achieving state-of-the-art performance on CIFAR-$10$ is at least $1.5$ million [4]. These neural networks have the ability to memorize pure noise – somehow, they are still able to generalize well when trained with meaningful data.

Since the publication of this observation [11, 4], the machine learning community has seen a flurry of activity to attempt to explain this phenomenon, both for classification and regression problems, in neural networks. The problem is challenging for three core reasons222This exposition is inspired by Suriya Gunasekar’s presentation at the Simons Institute, Summer $2019$.:

The optimization landscape for loss functions on neural networks is notoriously non-convex and complicated, and even proving convergence guarantees to some global minimum, as is observed in practice in the overparameterized regime [11, 4], is challenging. 2. 2.

In the overparameterized regime, there are multiple global minima corresponding to a fixed neural network architecture and loss function – which of these minima the optimization algorithm selects is not always clear. 3. 3.

Tight generalization bounds for the global minimum that is selected need to be obtained to show that overparameterization can help with generalization. This is particularly non-trivial to establish for deep, i.e. $\geq 3$-layer neural networks.

Promising progress has been made in all of these areas, which we recap only briefly below. Regarding the first point, while the optimization landscape for deep neural networks is non-convex and complicated, several independent recent works (an incomplete list is [12, 13, 14, 15, 16, 17, 18]) have shown that overparameterization can make it more attractive, in the sense that optimization algorithms like stochastic gradient descent (SGD) are more likely to actually converge to a global minimum. These interesting insights are mostly unrelated to the question of generalization, and should be viewed as a coincidental benefit of overparameterization.

Second, a line of recent work [19, 20, 21, 22] characterizes the inductive biases of commonly used optimization algorithms, thus providing insight into the identity of the global minimum that is selected. For the logistic loss of linear predictors on separable data, the global optimum that is found by SGD in overparameterized settings has been shown to be precisely the margin-maximizing support vector machine (SVM) [19], and extensions have been obtained for deep linear networks [21] and stochastic mirror descent [20]. More pertinent to linear regression, the behavior of SGD on the quadratic loss function [22] depends on the initialization point: the global convergence guarantee smoothly interpolates between a “kernel" regime (corresponding to the minimum $\ell_{2}$-norm interpolator) and a “structured" regime (corresponding to the minimum $\ell_{1}$-norm interpolator).

Finally, simple theoretical insights into which solutions (global minima) generalize well under what conditions, if any, remain elusive. The generalizing ability of solutions can vary for different problem instances: adaptive methods333for which, interestingly, the induced inductive bias is unknown. in optimization need not always improve generalization for specially constructed examples in overparameterized linear regression [23]. On the other hand, on a different set of examples [24], adaptive methods converge to better-generalizing solutions than SGD. Evidence suggests that norm-based complexity measures predict generalizing ability [11], and for neural networks, such complexity measures have been developed that do not depend on the width, but can depend on the depth [25, 26, 27, 28]. A classification-centric explanation for the possibility of overparameterization improving generalization is that SGD with logistic-style loss converges to the solution that maximizes training data margin. Then, increasing overparameterization increases model flexibility, allowing for solutions that increase the margin. This is a classical observation for AdaBoost [29] and is recently given as a justification for using overparameterization in neural networks [30]. However, margin does not always imply generalization [24]. An alternative, intriguing explanation [31] for this phenomenon does not consider margin, but instead connects the AdaBoost procedure to random forests, i.e. ensembles of randomly initialized decision trees, each of which interpolate the training data. An averaging and localizing-of-noise effect that is shown to be present both in the random forests ensemble, and the interpolating ensemble of AdaBoost, results in good generalization.

1.1.2 Kernels for interpolation of data

In the overparameterized regime, solutions that minimize classification/regression loss interpolate the training data. Another class of functions that have the ability to interpolate the training data are not explicitly overparameterized – they are non-parametric functions corresponding to particular reproducing kernel Hilbert spaces (RKHS). In fact, Belkin, Ma and Mandal [32] empirically recovered several of the overparameterization phenomena of deep learning in kernel classifiers that interpolate444In the paper, a subtle distinction is made between overfitting, which corresponds to close to zero classification loss, and interpolation, which corresponds to close to zero squared loss. the training data. They observed that these solutions generalize well, even in the presence of label noise; and moreover, regularization (either through explicit norm control or early stopping of SGD) yields only a marginal improvement555This was also observed by [11] with overparameterized neural networks.. Finally, they showed that the minimum (RKHS) norm of such interpolators in the presence of label noise cannot explain these properties; and generalizing ability likely depends on specific structure of the kernel. Other interpolators (e.g. based on local methods) have subsequently been analyzed for specific kernels [33, 34], but most relevant to the setting of regression is recent analysis of the test MSE of the minimum-RKHS-norm kernel interpolator [35]. It was shown here that this could be controlled under appropriate conditions on the eigenvalues of the kernel as well as the data matrix, and critical to these results is the dimension of the data growing with the number of samples [36]. Implicitly, all the analyses require successful interpolators to have a delicate balance between preserving the structure of the true function explaining the (noiseless) data and minimizing the harmful effect of regression noise in the data: thus, the properties of the chosen kernel are key. We will see that this tradeoff manifests very explicitly and clearly in high-dimensional linear regression.

1.1.3 High-dimensional linear regression

Most recently, a double-descent curve on the test error ($0-1$ loss and MSE) as a function of the number of parameters of several parametric models was observed on several common datasets by physicists [37] and machine learning researchers [38] respectively. In these experiments, the minimum $\ell_{2}$-norm interpolating solution is used, and several feature families, including kernel approximators [39], were considered. The experiments showed that the effect of increased parameterization on the test MSE is nuanced: at $d\sim n$, fitting noise has an adverse effect on the test MSE; but as the number of features $d$ becomes many times more than the number of data points $n$, the test MSE decays and approaches the test MSE of a corresponding kernel interpolator of minimum RKHS norm.

When the number of features is greater than the number of training data points, this corresponds to overparameterized linear regression (on lifted features). This is more commonly known as high-dimensional, or underdetermined linear regression. The classical underdetermined analysis [5] tells us that the data matrix could be ill-conditioned, especially at $d\sim n$: thus, underdetermined linear regression is usually carried out with some form of regularization. One of the earliest regularization procedures is Tikhonov ($\ell_{2}$) regularization, which can be used to improve the conditioning of the data matrix - elegant analysis of this for a random design, but in the overdetermined regime, is considered in [40]. In high dimensions, $\ell_{2}$-regularization unfortunately leads to a generic “bleeding" across features that causes a loss of the signal; this classical observation, and exceptions to this observation, are discussed in detail in Section 4. In fact, signal recovery in the high-dimensional regime, even in the absence of noise, is only possible when there is sparsity in the coefficients666This is precisely shown by information-theoretic lower bounds [41, 42] and justifies the principle to “bet on sparsity” [1].. The literature on sparse signal recovery is rich, with algorithms that are motivated by iterative perspectives from signal processing (orthogonal matching pursuit (OMP) [43] and stagewise OMP [44]) and convex relaxation (basis pursuit/Lasso [45]). Unlike $\ell_{2}$-regularized solutions, which tend to spread energy out across multiple features, these procedures are parsimonious and select features in a data-dependent manner. Traditionally, these solutions are applied on denoised data – or, equivalently, with regularization777To be precise, they are analyzed for sufficiently large values of the regularizer $\lambda$. An analysis of the equivalent interpolators corresponds to an analysis tending $\lambda\to 0$. to explicitly avoid fitting noise. We investigate the performance of solutions that interpolate both signal and noise in the high-dimensional, sparse linear model in Section 5.

1.1.4 Concurrent work in high-dimensional linear regression

Since the publication888In fact, as pointed out in [8], the initial discussions around this area began at the Simons Institute $2017$ program on foundations of machine learning. of the double descent experiments [37, 38] at the end of the year $2018$, there has been extremely active interest in understanding the generalization abilities of interpolating solutions in linear regression. An earlier edition of our work was presented at Information Theory and Applications, February $2019$ and subsequently accepted to IEEE International Symposium on Information Theory, July $2019$. Several elegant and interesting papers [6, 8, 7, 9, 10, 46] have appeared around this time. All of these center around the analysis of the $\ell_{2}$-minimizing interpolator. We discuss [6, 8, 7] in substantial detail in Section 4, and provide a succinct description of each of these paper’s contributions here:

Belkin, Hsu and Xu [7] consider an effectively misspecified setting for Gaussian and Fourier features and recover the double descent phenomenon for the $\ell_{2}$-minimizing interpolator in cases where the data can be generated as a linear combination of a huge number of features – therefore, adding more features into the model in a “prescient" manner can have an approximation-theoretic benefit. On the other hand, double descent no longer happens under this model when the features are randomly selected. 2. 2.

Hastie, Tibshirani, Rosset and Montanari [6] analyze the asymptotic risk of the $\ell_{2}$-minimizing interpolator for independent and whitened features with a bounded fourth moment as a function of the “overparameterizing" factor: $\frac{\text{number of features}}{\text{number of samples}}$. They asymptotically prove several phenomena corresponding to the non-asymptotic results in our paper – harmless noise fitting with increased overparameterizing factor (our Corollaries 1 and 2), what we call (in our discussion in Section 4) signal “bleed" and the resulting convergence to the null risk. Additionally, they study the misspecified model asymptotically and show that double descent can occur with sufficiently improved approximability, and replicate the phenomena in a specific non-linear model. 3. 3.

Bartlett, Long, Lugosi and Tsigler [8] sharply upper and lower bound the (non-asymptotic) generalization error of the $\ell_{2}$-minimizing interpolator for Gaussian features (which are not whitened/independent in general). They characterize necessary and sufficient conditions for the $\ell_{2}$-minimizing interpolator to avoid what we call signal “bleed" and noise overfitting in terms of functionals of the spectrum of the Gaussian covariance matrix. Initial discussions with Peter Bartlett as well as the paper’s random-matrix-theory-based analysis significantly inspired our subsequent Fourier-theoretic interpretation that we provide in Section 4. 4. 4.

Bibas, Fogel and Feder [9] study online/universal linear regression in the overparameterized regime and show that when the high-dimensional training data is mostly along certain pre-determined directions that favor the true parameter vector, good generalization is possible with both $\ell_{2}$-regularized and minimum $\ell_{2}$-norm interpolators used for prediction. 5. 5.

Mitra [10] also studies the asymptotically overparameterized regime and obtains precise analytical risk curves for both the $\ell_{2}$-minimizing interpolator and the $\ell_{1}$-minimizing interpolator, with special focus on the overfitting peak at the “interpolation threshold", i.e. $d\sim n$.

Subsequent to all of this work, Mei and Montanari [46] very recently proved that Tikhonov regularization in the overparameterized regime can give rise to the double-descent behavior when random Fourier features are used, in both the noiseless and noisy cases. They also recover this behavior with $\ell_{2}$-minimizing interpolation, but in the low-noise (high SNR) limit. For a more complete understanding of interpolation in overparameterized linear regression, we recommend to the reader all of these papers along with ours.

2 Problem Setting

Throughout, we consider data that is actually generated from the high-dimensional/overparameterized linear model999The more general misspecified case is considered in concurrent and subsequent work [6, 46].. We consider covariate-response pairs $(\mathbf{X}_{i},Y_{i}\in\mathbb{R}^{p}\times\mathbb{R})_{i=1}^{n}$ and generative model $Y=\langle\mathbf{a}(\mathbf{X}),\,\bm{\alpha}^{*}\rangle+W$ for feature vector $\mathbf{a}(\mathbf{X})\in\mathbb{R}^{d}$ and Gaussian noise $W\sim\mathcal{N}(0,\sigma^{2})$ that is independent of $X$. We generically assume that the covariates $\{\mathbf{X}_{i}\}_{i=1}^{n}$ are iid random samples, but will also consider regularly spaced data on bounded domains for some special cases. The signal $\bm{\alpha}^{*}$ is unknown apriori to an estimator. We also assume a distribution on $\mathbf{X}\in\mathbb{R}^{p}$, which induces a distribution on the $d$-dimensional feature vector $\mathbf{a}(\mathbf{X})$. Let $\bm{\Sigma}=\mathbb{E}[\mathbf{a}(\mathbf{X})\mathbf{a}(\mathbf{X})^{\top}]$ denote the covariance matrix of the feature vector under this induced distribution. We assume that $\bm{\Sigma}$ is invertible; therefore it is positive definite and its square-root-inverse $\bm{\Sigma}^{-1/2}$ exists.

We define shorthand notation for the training data: let

[TABLE]

denote the data (feature) matrix, and let $\mathbf{Y}_{\mathsf{train}},\mathbf{W}_{\mathsf{train}}\in\mathbb{R}^{n}$ denote the output and noise vectors respectively.

We will primarily consider the overparameterized, or high-dimensional regime, i.e. where $d>n$. We are interested in solutions $\bm{\alpha}$ that satisfy the following feasibility condition for interpolation:

[TABLE]

We assume that $\text{rank}(\mathbf{A}_{\mathsf{train}})=n$, so the set $\{\bm{\alpha}\in\mathbb{R}^{d}:\mathbf{A}_{\mathsf{train}}\bm{\alpha}=\mathbf{Y}_{\mathsf{train}}\}$ is non-empty in $\mathbb{R}^{d}$.

For any solution $\widehat{\bm{\alpha}}\in\mathbb{R}^{d}$, we define the generalization error as test MSE101010This constitutes all quadratic loss functions on error in estimation of the signal $\bm{\alpha}^{*}$, and is the standard error metric for regression. below.

Definition 1**.**

The expected test mean-squared-error (MSE) minus irreducible noise error of any estimator $\widehat{\bm{\alpha}}((\mathbf{X}_{i},Y_{i})_{i=1}^{n})$ is given by

[TABLE]

where the expectation is taken only over the joint distribution on the fresh test sample $(\mathbf{X},Y)$, and we subtract off the irreducible noise error $\mathbb{E}[W^{2}]=\sigma^{2}$.

We have chosen the convention to subtract off the unavoidable error arising from noise, $\sigma^{2}$, as is standard. From now on, we will denote this quantity to be the test MSE as shorthand.

3 The fundamental price of interpolation

Before analyzing particular interpolating solutions, we want to understand whether interpolation can ever lead to a desirable guarantee on the test MSE $\mathcal{E}_{\mathsf{test}}$. To do this, we characterize the fundamental price that any interpolating solution needs to pay in test MSE. The constraint in Equation (1) is sufficiently restrictive to not allow trivial solutions of the form $\bm{\alpha}=\bm{\alpha}^{*}$ — so this is a surprisingly well-posed problem. In fact, we can easily define the ideal interpolator below.

Definition 2**.**

The ideal interpolator $\widehat{\bm{\alpha}}_{\mathsf{ideal}}$ is defined as:

[TABLE]

We also denote the test MSE of the ideal interpolator, which we henceforth call the ideal test MSE, as $\mathcal{E}_{\mathsf{test}}^{*}$. This is, by definition, a lower bound on the test MSE of any interpolator.

The following result exactly characterizes the ideal interpolator and the ideal test MSE.

Theorem 1**.**

For any joint distribution on $(\mathbf{X},Y)$ and realization of training data matrix $\mathbf{A}_{\mathsf{train}}$ and noise vector $\mathbf{W}_{\mathsf{train}}$, the lowest possible test MSE any interpolating solution can incur is bounded below as $\mathcal{E}_{\mathsf{test}}\geq\mathcal{E}_{\mathsf{test}}^{*}$, where

[TABLE]

Here, $\mathbf{B}_{\mathsf{train}}:=\mathbf{A}_{\mathsf{train}}\bm{\Sigma}^{-1/2}$ is the whitened training data matrix.

The proof of Theorem 1 is outlined in Section 3.3. Theorem 1 provides an explicit expression for a fundamental limit on the generalization ability of interpolation. Thus, we can easily evaluate it (numerically) when the training data matrix $\mathbf{A}_{\mathsf{train}}$ is generated by a number of choices for feature families $\mathbf{a}(\mathbf{X})$. These choices are listed below as examples.

Example 1** (Gaussian features).**

The Gaussian features on $d$-dimensional data comprise of $\mathbf{a}(\mathbf{X}):=\mathbf{X}\sim\mathcal{N}(\mathbf{0},\bm{\Sigma})$, where $\bm{\Sigma}\in\mathbb{R}^{d\times d}$ and $\bm{\Sigma}\succ 0$. A special case is iid Gaussian features, i.e. $\bm{\Sigma}=\mathbf{I}_{d}$.

Example 2** (Fourier features).**

Let $i:=\sqrt{-1}$ denote the imaginary number. For one-dimensional data $X\in[0,1]$, we can write the $d$-dimensional Fourier features in their complex form as

[TABLE]

*This is clearly an orthonormal feature family in the sense that

$\mathbb{E}_{X\sim\text{Unif}[0,1]}\left[\mathbf{a}(X)_{j}\mathbf{a}(X)_{k}^{*}\right]=\delta_{j,k}$, where $\delta_{j,k}$ denotes the Kronecker delta and $(\cdot)^{*}$ denotes the complex conjugate. When evaluating interpolating solutions for Fourier features, we will consider one of two models for the data $\{X_{i}\}_{i=1}^{n}$:*

$n$-regularly spaced training data points, i.e. $x_{i}=\frac{(i-1)}{n}\text{ for all }i\in[n]$, which we consider empirically and theoretically in Section 4. 2. 2.

$n$-random training data points, i.e. $X_{i}\text{ i.i.d }\sim\text{Unif}[0,1]$, which we evaluate only empirically.

Example 3** (Legendre/Vandermonde polynomial features).**

For one-dimensional data $X\in[-1,1]$, we can write the $d$-dimensional Vandermonde features as

[TABLE]

We can also uniquely define their orthonormalization with respect to the uniform measure on $[-1,1]$. In other words, we define the $d$-dimensional Legendre features as polynomials

[TABLE]

where $\text{deg}(p_{j}(X))=j$ for every $j\geq 0$, and $\{p_{j}(X)\}_{j\geq 0}$ is uniquely defined such that $\mathbb{E}_{X\sim\text{Unif}[-1,1]}\left[p_{j}(X)p_{k}(X)\right]\allowbreak=\delta_{j,k}$, i.e. the Legendre polynomials form the orthonormal basis with respect to the uniform measure on $[-1,1]$. When evaluating interpolating solutions for both these polynomial features, we consider one of two models for the training data $\{X_{i}\}_{i=1}^{n}$:

$n$-regularly spaced training data points, i.e. $x_{i}=-1+\frac{2(i-1)}{n}\text{ for all }i\in[n]$. 2. 2.

$n$-random training data points, i.e. $X_{i}\text{ i.i.d }\sim\text{Unif}[-1,1]$.

Figure 7 evaluates the quantity $\mathcal{E}_{\mathsf{test}}^{*}$ as a function of $d$ for iid Gaussian features. For $d<n$, we always evaluate the test MSE of the unique least-squares solution as Equation (1) is no longer feasible. We observe a spike at $d=n$, and a decay in the generalization error as $d>>n$, implying that potentially harmful effects of noise can be mitigated for these feature families. These properties also manifest in Figure 1 for the orthonormalized Legendre polynomial features, but not for the Vandermonde polynomial features — illustrating the importance of the whitening step.

3.1 Converse: Lower bound on test MSE for interpolating solutions

As mentioned in Definition 2, the quantity $\mathcal{E}_{\mathsf{test}}^{*}$ represents the minimum test error that any interpolating solution satisfying Equation (1) must incur. Noteworthy from the expression in Equation (2) is that the test MSE naturally depends on the singular value spectrum of the random matrix $\mathbf{B}_{\mathsf{train}}^{\top}$ that has columns $\{\mathbf{b}(\mathbf{X}_{i})\}_{i=1}^{n}$ corresponding to whitened features, i.e. $\mathbb{E}\left[\mathbf{b}(\mathbf{X}_{i})\mathbf{b}(\mathbf{X}_{i})^{\top}\right]=\mathbf{I}_{d}$ for every $i\in[n]$. Because of this whitening, we might expect that the matrix $\mathbf{B}_{\mathsf{train}}^{\top}$ becomes better and better conditioned and the error arising solely from fitting noise reduces. To understand the fundamental price of interpolation, we make one of three assumptions on columns of $\mathbf{B}_{\mathsf{train}}^{\top}$, denoted by $\{\mathbf{b}(\mathbf{X}_{i}):=\bm{\Sigma}^{-1/2}\mathbf{a}(\mathbf{X}_{i})\}_{i=1}^{n}$. Recall that the random feature vectors $\mathbf{b}(\mathbf{X}_{i})$ are always iid, because we have assumed iid samples $\{\mathbf{X}_{i}\}_{i=1}^{n}$. We state our candidate assumptions in increasing order of strength.

Assumption 1**.**

*[iid heavy-tailed feature vectors] The random feature vectors $\mathbf{b}(\mathbf{X}_{i})$ are bounded almost surely, i.e. $\|\mathbf{b}(\mathbf{X}_{i})\|_{2}\leq\sqrt{d}$ almost surely. Note that this assumption is satisfied by discrete Fourier features and random Fourier features. *

Assumption 2**.**

[iid feature vectors with sub-Gaussianity] The whitened feature vectors $\mathbf{b}(\mathbf{X}_{i})$ are sub-Gaussian with parameter at most $K>0$. A special case of this includes the case of independent entries: in this case, the random variables $B_{ij}$ are independent and sub-Gaussian, all with parameter at most $K>0$. We reproduce the definition of sub-Gaussianity and sub-Gaussian parameter for both random vector and random variable in Definition 6 in Appendix A.1.

Assumption 3**.**

[iid Gaussian entries] The entries of the whitened feature matrix $\mathbf{B}_{\mathsf{train}}$ are iid Gaussian, i.e. $B_{ij}\text{ i.i.d }\sim\mathcal{N}(0,1)$. Note that this exactly describes all cases where the original data matrix $\mathbf{A}_{\mathsf{train}}$ has iid Gaussian row vectors, i.e. $\mathbf{a}(\mathbf{X}_{i})\sim\mathcal{N}(\mathbf{0},\Sigma)$.

Notice that the Gaussian Assumption 3 constitutes a special case of sub-Gaussianity of rows (Assumption 2). Independence of elements of the feature vector, even when the features are whitened, is impossible when lower-dimensional data is lifted into high-dimensional features, i.e. the problem is one of lifted linear regression. It is in view of this that we have included consideration of the far weaker assumptions of sub-Gaussianity of random feature vectors (Assumption 2) and even heavy-tailed features (Assumption 1). We will see that the strength of the conclusions we can make is accordingly lower for these more general cases. However, for random feature vectors satisfying any of the above assumptions (which, together, constitute very mild conditions), we can always characterize the fundamental price of interpolation by lower bounding the ideal MSE.

Corollary 1**.**

For any $0<\delta<1$, the fundamental price of any interpolating solution is at least:

[TABLE]

with probability greater than or equal to $(1-\frac{1}{n}-e^{-\delta^{2}n/8})$ for any feature family for which the random whitened feature matrix $\mathbf{B}_{\mathsf{train}}$ satisfies the heavy-tailed Assumption 1. Here, $C>0$ is some positive constant independent of the choice of feature family. 2. 2.

[TABLE]

*with probability greater than or equal to $(1-e^{-c_{K}n}-e^{-\delta^{2}n/8})$ for any feature family for which the random whitened feature matrix $\mathbf{B}_{\mathsf{train}}$ satisfies Assumption 2 of iid sub-Gaussian rows. (This includes the special case in which $\mathbf{B}_{\mathsf{train}}$ has independent sub-Gaussian entries.) Here, $C_{K},c_{K}>0$ are positive constants that depend on the upper bound on the sub-Gaussian parameter, $K>0$. * 3. 3.

[TABLE]

with probability greater than or equal to $(1-e^{-n/2}-e^{-n\delta^{2}/8})$ for any Gaussian feature family, i.e. any feature family satisfying the Gaussian Assumption 3.

Corollary 1 characterizes the fundamental price of any interpolating solution as $\widetilde{\mathcal{O}}\left(\frac{n}{d}\right)$ at a significant level of generality111111Note that in the case of data satisfying Assumption 1, the $\widetilde{\mathcal{O}}(\cdot)$ omits the $(\ln n)$ factor in the denominator arising from heavy tailed-ness.. It tells us that extreme overparameterization is essential for harmlessness of interpolation of noise. To see this, consider how the number of features $d$ could scale as a function of the number of samples $n$. Say that we grew $d=\gamma n$ for some constant $\gamma>1$. Then, the lower bound on test MSE (minus the irreducible error $\sigma^{2}$ arising from the prospect of noise in the test points as well) scales as $\omega\left(\frac{\sigma^{2}n}{d}\right)=\omega\left(\frac{\sigma^{2}}{\gamma}\right)$, which asymptotes to a constant as $n\to\infty$. This tells us that the level of overparameterization necessarily needs to grow faster than the number of samples for harmless interpolation to even be possible. For example, this could happen at a polynomial rate (e.g. $d=n^{q}$ for some $q>1$) or even an exponential rate (e.g. $d=e^{\lambda n}$ for some $\lambda>0$).

3.2 The possibility of harmless interpolation

Corollary 1 only provides a lower bound on the ideal MSE, not an upper bound. Thus, it does not tell us whether harmless interpolation is ever actually possible. This turns out to be a more delicate question in general, and is difficult to characterize under the weaker Assumptions 1 and 2 in their full generality (for a detailed discussion of why this is the case, see Appendix A.1). However, for two special cases: a) Gaussian features (Assumption 3), and b) independent sub-Gaussian feature vectors (Assumption 2 for independent entries): we can show that harmless interpolation is always possible, with an upper bound on the ideal MSE that matches the lower bound provided in Corollary 1. We state this result below.

Corollary 2**.**

For any $0<\delta<1$ and:

Random whitened feature matrix $\mathbf{B}_{\mathsf{train}}$ satisfying the Gaussian Assumption 3, the fundamental price of interpolation is at most

[TABLE]

with probability greater than or equal to $(1-e^{-n/2}-e^{-n\delta^{2}/8})$. 2. 2.

Random matrix feature matrix $\mathbf{B}_{\mathsf{train}}$ satisfying independent, sub-Gaussian entries with unit variance (special case of Assumption 2), the fundamental price of interpolation is at most

[TABLE]

with probability greater than or equal to $\left(1-\left(\frac{1}{2}\right)^{d-n+1}-c_{K}^{d}-e^{-n\delta^{2}/8}\right)$, where constants $c_{K}\in(0,1)$ and $C_{K}>0$ only depend on the upper bound on the sub-Gaussian parameter $K$.

Thus, when the features given to us in the form of training data matrix $\mathbf{A}_{\mathsf{train}}$ are normally distributed (in general unwhitened), Corollary 2 helps explain the harmless interpolations that we are seeing empirically. Together with Corollary 1, this is a sharp, up to constants121212The concurrent work of Hastie, Tibshirani, Rosset and Montanari [6] exactly characterizes the effect of noise-fitting for feature vectors with independent entries (Lemma $3$ in their paper) satisfying a bounded fourth moment assumption in the asymptotic regime where $d,n\to\infty,\lim\frac{d}{n}=\gamma$. Their scaling in $\gamma$ is equivalent to ours in the constant high-dimensional regime $d=\gamma n$. Their result is a consequence of being able to exactly characterize the asymptotic distribution of the spectrum of the infinite-dimensional random matrix $\mathbf{B}_{\mathsf{train}}\mathbf{B}_{\mathsf{train}}^{\top}$ by applying a generalized Marchenko-Pastur law. We use non-asymptotic random matrix theory for our results and thus sharply characterize the dependence on $n$ and $d$ in our rates, but not the dependence on constant factors. , characterization of the fundamental price of any interpolating solution when using a Gaussian feature family131313For the Gaussian case, this scaling can also be derived from the work of Bartlett, Long, Lugosi and Tsigler [8], which analyes the generalization error of the $\ell_{2}$-minimizing interpolation of unwhitened Gaussian features to Gaussian noise more generally. As we see from the second part of Corollary 2, this scaling holds more generally for whitened, iid sub-Gaussian features. . To make this notion of harmlessness more concrete, it is useful to consider two ultra-overparameterized settings:

$d=n^{q}$ for some parameter $q>1$, which constitutes a polynomially high-dimensional regime. Here, the ideal test MSE is upper bounded by $\mathcal{O}\left(\frac{\sigma^{2}n}{n^{q}}\right)=\omega\left(\frac{\sigma^{2}}{n^{1-q}}\right)$, which goes to [math] as $n\to\infty$ at a polynomial rate. 2. 2.

$d=e^{\lambda n}$ for some parameter $\lambda>0$, which constitutes an exponentially high-dimensional regime. Here, ideal test MSE is at most $\mathcal{O}\left(\sigma^{2}ne^{-\lambda n}\right)$ which goes to [math] as $n\to\infty$ at an exponentially decaying rate.

Thus, there always exists an interpolating solution that fits noise in such a manner that the effect of fitting this noise on test error decays to [math] as the number of features goes to infinity. Of course, Corollary 2 is not particularly meaningful for $d\sim n$, i.e. near the interpolation threshold: for a detailed discussion of this regime, see Section B.

We defer the proofs of Corollary 1 and Corollary 2 for the more general cases of sub-Gaussian and heavy-tailed feature vectors (Assumptions 1 and 2) to Appendix A.1. We here prove Theorem 1 and Corollaries 1 and 2 for the iid Gaussian case (Assumption 3).

3.3 Proof of Theorem 1, Corollaries 1 and 2

To prove Theorem 1, we first get an exact expression for the ideal interpolator as defined in Definition 2. A simple calculation gives us

[TABLE]

Thus, we can equivalently characterize the ideal interpolating solution $\widehat{\bm{\alpha}}_{\mathsf{ideal}}:={\arg\min}\text{ }\mathcal{E}_{\mathsf{test}}(\widehat{\bm{\alpha}})$ as:

[TABLE]

Observe that Equation (1) can be rewritten as

[TABLE]

Then we have a closed form expression for the minimum norm solution, denoting $\mathbf{B}_{\mathsf{train}}=\mathbf{A}_{\mathsf{train}}\bm{\Sigma}^{-1/2}$:

[TABLE]

(note that $\mathbf{B}_{\mathsf{train}}^{\dagger}:=\mathbf{B}_{\mathsf{train}}^{\top}(\mathbf{B}_{\mathsf{train}}\mathbf{B}_{\mathsf{train}}^{\top})^{-1}$ is just the right Moore-Penrose pseudoinverse of $\mathbf{B}_{\mathsf{train}}$).

Substituting this expression into the test MSE calculation:

[TABLE]

Plugging this into the expression of $\mathcal{E}_{\mathsf{test}}^{*}$ gives us Equation (2), thus proving Theorem 1. ∎

To prove Corollaries 1 and 2, we lower bound and upper bound the test MSE of the ideal interpolator in Equation (2). We use matrix concentration theory [47] to do this for the Gaussian case (Assumption 3). We first state the following concentration result on the non-zero singular values of random matrix $\mathbf{B}$ with entries $B_{ij}\text{ i.i.d }\sim\mathcal{N}(0,1)$ as is from Vershynin’s book [47]. The original argument is contained in classical work by Davidson and Szarek [48].

Lemma 1** (Theorem $5.32$,  [47]).**

Let $\sigma_{min}(\cdot)$ and $\sigma_{max}(\cdot)$ denote the minimum and maximum (non-degenerate) singular values of a matrix. For random matrix $\mathbf{B}\in\mathbb{R}^{n\times d}$ having entries $B_{ij}\text{ i.i.d. }\sim\mathcal{N}(0,1)$, we have

[TABLE]

with probability at least $(1-e^{-t^{2}/2})$.

We first use this lemma to prove Corollary 1. We lower bound Equation (2) as

[TABLE]

Now, we apply the upper bound on the maximum singular value of $\mathbf{B}_{\mathsf{train}}^{\top}$ as stated in Lemma 1, substituting $t:=\sqrt{n}$ to get

[TABLE]

with probability at least $(1-e^{-n/2})$. Further, we have $\|\mathbf{W}_{\mathsf{train}}\|_{2}^{2}=\sum_{i=1}^{n}W_{i}^{2}$ and the lower tail bound on chi-squared random variables [49, Chapter $2$] gives us

[TABLE]

with probability greater than or equal to $(1-e^{-n\delta^{2}/8})$ for any $0<\delta<1$. When Equations (9) and (10) both hold, we get the statement of Corollary 1. ∎

Now, we prove Corollary 2. Denoting $\mathbf{B}_{\mathsf{train}}^{T}=\begin{bmatrix}\mathbf{b}_{1}&\ldots&\mathbf{b}_{n}\end{bmatrix}$, we can upper bound Equation (2) as

[TABLE]

Observe that the matrix $\mathbf{B}_{\mathsf{train}}$ has entries $B_{ij}\sim\mathcal{N}(0,1)$ due to the whitening. Thus, we can apply the lower bound on the minimum singular value of $\mathbf{B}_{\mathsf{train}}^{\top}$ as stated in Lemma 1, substituting $t:=\sqrt{n}$ to get

[TABLE]

with probability at least $(1-e^{-n/2})$.

Further, we have $\|\mathbf{W}_{\mathsf{train}}\|_{2}^{2}=\sum_{i=1}^{n}W_{i}^{2}$ and the corresponding upper tail bound on chi-squared random variables141414One could have also used the Hanson-Wright inequality [50] to get slightly more precise constants, but the tight concentration of the singular values of the random matrix $\mathbf{B}_{\mathsf{train}}$ implies that only constant factors would be improved. gives us

[TABLE]

with probability greater than or equal to $(1-e^{-n\delta^{2}/8})$ for any $\delta>0$. When Equations (11) and (12) both hold, we get the upper inequality in Equation (6), thus proving Corollary 2.

Using the union bound on the probability of non-event of Equations (9), (10), (11) or (12), gives us the statements of Corollaries 2 as well as 1 with probability greater than or equal to $(1-2e^{-n/2}-2e^{-\delta^{2}n/8})$. This provides a characterization of the fundamental price of interpolation for the Gaussian case. ∎

The proof of Theorem 1 tells us that the ideal interpolator fits what can be thought of as a residual signal, $(\bm{\alpha}-\bm{\alpha}^{*})$, to noise. This is well and truly an ideal interpolator, as it requires knowledge of the true signal $\bm{\alpha}^{*}$ to implement – and thus it should only be thought of as primarily a fundamental limit on generalization. To understand whether we could achieve this limit, we turn to practically used interpolating schemes, beginning with the $\ell_{2}$-minimizing interpolator which has seen the most concurrent (empirical and theoretical) attention.

4 The minimum-$\ell_{2}$-norm interpolator through the Fourier lens

For overparameterized linear models, it is natural to consider the minimum-$\ell_{2}$-norm solution that interpolates the data, i.e. satisfies Equation (1). After all, the proof of Theorem 1 shows that the least harmful way in which to fit noise is to find the minimum-$\ell_{2}$-norm interpolator of (effectively) pure noise using appropriately whitened features. Furthermore, the minimum-$\ell_{2}$-norm interpolator is explicitly characterizable as a linear matrix operator on the output vector $\mathbf{Y}_{\mathsf{train}}$, and can be easily computed as well. For example, gradient descent on the overparameterized linear regression problem with appropriate initialization converges to the minimum-$\ell_{2}$-norm interpolator [24, 22].

As mentioned in Section 1.1, a couple of excellent papers [6, 8] that center around comprehensive analyses of the $\ell_{2}$-minimizing interpolator were published concurrently151515Through personal communication with Peter Bartlett in December 2018/January 2019, we became aware of the gist of the investigations underlying [8] when we first presented our work. In the earlier version of this work presented at IEEE ISIT 2019, we cite this communication accordingly.. When it comes to whitened or Gaussian features, any results discussed in this section are all reflected, more or less, across [6] and [8] respectively – and we recommend reading these papers for the details. The results, like our proofs of Theorem 1, Corollary 1 and 2, use fundamental advances in asymptotic and non-asymptotic random matrix theory. Here, we provide a brief alternative exposition of the main ideas through a Fourier-theoretic lens on regularly spaced data (Example 2). Our aim for doing this is two-fold: one, simply to complement these papers; and two, to provide a signal processing oriented perspective on salient properties of the minimum-$\ell_{2}$-norm interpolator.

From Corollaries 1 and 2, we saw that the fundamental price of interpolation of noise on test MSE scales as $\Theta\left(\frac{n}{d}\sigma^{2}\right)$, where $n$ is the number of training samples, $d$ is the number of features, and $\sigma^{2}$ is the noise variance. Looking at Figure 2, we see that when $d$ is large, this is also the scaling that is achieved by the case of regularly spaced data points with Fourier features. This case, as an easy-to-understand paradigmatic example, provides a clear lens into understanding what is happening conceptually *for $\ell_{2}$-minimizing solutions161616What we will use is the exact presence of exactly aliased higher-frequency features corresponding to any low-frequency features. This will give extremely clean expressions for the adverse effect that $\ell_{2}$-minimization has on recovering a signal with low-frequency components. It is not as simple to understand the effect of sparsity-seeking methods, which require more generic properties like restricted isometry to work (although such properties have also been established for such sub-sampled Fourier matrices under nonuniform random sampling [51]). *. It is first useful see how the minimum-$\ell_{2}$-norm interpolator actually behaves for two contrasting examples.

Example 4** (Standard Gaussian features, $k=500$-sparse signal).**

We consider $d$-dimensional iid standard Gaussian features, i.e. Example 1 with $\bm{\Sigma}=\mathbf{I}_{d}$. In other words, the features $\{a(\mathbf{X})_{j}\}_{j=1}^{d}$ are iid and distributed as $\mathcal{N}(0,1)$. Let the first $500$ entries of the true signal $\bm{\alpha}^{*}$ be non-zero and the rest be zero, i.e. $\mathsf{supp}(\bm{\alpha}^{*})=[500]$. We take $n=5000$ measurements, each of which is corrupted by Gaussian noise of variance $\sigma^{2}=0.01$.

Figure 7 shows the test MSE of the minimum-$\ell_{2}$-norm interpolator on Example 4 as a function of the number of features $d$ for various choices of noise variance $\sigma^{2}$. We immediately notice that the test MSE is converging to the same level as the test error from simply using a hypothetical $\bm{\alpha}=\mathbf{0}$. This property of generalization error degenerating to that of simply predicting [math] was also pointed out in [6] where this level was called out as the “null risk.”

Example 5** (Standard Gaussian features plus a constant, constant signal).**

In this example, we consider $d$-dimensional iid Gaussian features with unit mean and variance equal to $0.01$. More precisely, the features $\{a(\mathbf{X})_{j}\}_{j=1}^{d}$ are iid and distributed as $\mathcal{N}(1,0.01)$. We also assume the generative model for the training data:

[TABLE]

where as before, $W\sim\mathcal{N}(0,\sigma^{2})$ is the observation noise in the training data, and we pick $\sigma^{2}=0.01$. Note that in this example the true “signal" is the constant $1$, which is not exactly expressible as a linear combination of the $d$ Gaussian features. We take $n=10$ noisy measurements of this signal.

Example 5 is directly inspired by feature models171717A more general version of this example would replace the constant $1$ in the means of the features by the relevant realizations of an underlying latent Gaussian random feature vector with the true signal being that latent Gaussian feature vector. The qualitative behavior of double-descent will be retained, and we defer a formal discussion of this to a forthcoming paper. in several recent papers [7, 52, 8], as well as having philosophical connections to other recent papers [46, 6]. Figure 3 shows the performance of the minimum-$\ell_{2}$-norm interpolator on Example 5 as we increase the number of features. Here, we clearly see the double-descent behavior as the test MSE of the minimum-$\ell_{2}$-norm interpolator decreases with increased overparameterization. Note from Equation (13) that the true signal is not exactly representable by a linear combination of the random Gaussian features, and in fact Example 5 is an instance of the linear model being misspecified for any (finite) number of features $d$. Improved approximability from adding more features to the model partially explains the double descent behavior, but it is not the main reason. It turns out that we would still see the double-descent behavior with the minimum-$\ell_{2}$-norm interpolator if we added another feature that was always the constant $1$ (in which case the example would fall under the $(1,\sigma^{2})$-sparse linear model; see Definition 3.)

The minimum-$\ell_{2}$-norm interpolator generalizes well for Example 5, showing double descent behavior – but extremely poorly for Example 4. Why does one case fail while the other case works? Bartlett, Long, Lugosi and Tsigler [8] give an account of what is happening directly using the language of random matrix theory. Their paper defines a distinct pair of “effective ranks" using the spectrum of the covariance matrix $\bm{\Sigma}$ to state their necessary and sufficient conditions for this interpolator generalizing well. Classic core concepts in signal processing provide an alternative lens to understand these conditions. To use the Fourier-theoretic lens, we will naturally map the number of regularly spaced training samples $n$ to what is called the sampling rate in signal processing, and the number of features $d$ to what is called the bandwidth in signal processing. What we call overparameterization is what is called potential undersampling in the signal processing literature.

4.1 Aliasing — the core issue in overparameterized models

We have mapped overparameterization to undersampling of a true signal. The fundamental issue with undersampling of a signal is one of identifiability: infinitely many solutions, each of which correspond to different signal functions, all happen to agree with each other on the $n$ regularly spaced data points. These different signal functions, of course, disagree everywhere else on the function domain, so the true signal function is not truly reconstructed by most of them. This results in increased test MSE when such an incorrect function is used for prediction. Such functions that are different, but agree on the sampled points, are commonly called aliases of each other in signal processing language. Exact aliases naturally appear when the features are Fourier, as we see in the below example.

Example 6** (Fourier features with a constant signal).**

Denote $i:=\sqrt{-1}$ as the imaginary number. Consider the Fourier features as defined in complex form in Example 2 and regularly spaced input on the interval $[0,1)$, i.e. $x_{j}=\frac{j-1}{n}$ for all $j\in[n]$.

Suppose the true signal is equal to $1$ everywhere and the sampling model in the absence of noise is

[TABLE]

The estimator has to interpolate this data with some linear combination of Fourier features181818Why are we using complex features for our example instead of the real sines and cosines? Just because keeping track of which feature is an alias of which other feature is less notationally heavy for the complex case. The essential behavior would be identical if we just considered sines and cosines. $f_{k}(x)=e^{i2\pi kx}$ for $k=0\ldots d$.

A trivial signal function that agrees with Equation (14) at all the data points is the first (constant) Fourier feature: $f_{0}(x)=e^{i2\pi(0)x}=1$. It is, however, not the only one. The complex feature $f_{n}(x)=e^{i2\pi(n)x}$ will agree with $f_{0}(x)=1$ on all the regularly spaced points $\{x_{j}\}_{j=1}^{n}$ by the cyclic property of complex Fourier features (i.e. we have $e^{i2\pi(n)\frac{j-1}{n}}=e^{i2\pi(j-1)}=1=f_{0}(x)$). This is similarly true for features $f_{ln}(x)$ for all $l=1,2,\ldots,\frac{d}{n}-1$, and we thus have $\frac{d}{n}-1$ exact aliases191919These aliases are essentially higher frequency sinusoids that look the same as the low frequency one when regularly sampled at the rate $n$. This is the classic “movie of a fan under a strobelight” visualization where a fan looks like it is stationary instead of moving at a fast speed! of the true signal function $f_{0}(x)$ on the regularly spaced data points.

The above property is not unique to the constant function $f_{0}(x)$: for any true signal function that contains the complex sinusoid of frequency $k^{*}\in[n]$, i.e. $f_{k^{*}}(x)=e^{i2\pi(k^{*})x}$, the one-complete-cycle signal function $f_{k^{*}+n}=e^{i2\pi(k^{*}+n)x}$ again agrees on the $n$ regularly spaced data points, and for this signal function we again have the $\frac{d}{n}$ exact aliases $f_{k^{*}+ln}$ for all $l=1,2,\ldots,\frac{d}{n}-1$.

The presence of these aliases will naturally affect signal reconstruction. Before discussing this issue, we show an advantage in having aliases: they naturally absorb the noise that can harm generalization. Critically, as defined in Example 2, the Fourier features are orthonormal to each other in complex function space (where the integral that defines the complex inner product is taken with respect to the uniform measure on $[0,1)$). If we used only the first $n$ Fourier features (i.e. $f_{k}(x)$ upto $k=n-1$) to fit an $n$-dimensional pure noise vector (as described by $\{Y_{j}=W_{j}\}_{j=1}^{n}$), the coefficients of the $n$-dimensional fit, i.e. $\{\widehat{\alpha}_{k}\}_{k=1}^{n}$, would directly correspond to the appropriate discrete Fourier transform (DFT)202020The convention for defining the DFT depends on the chosen normalization. The symmetric/unitary DFT can be viewed as choosing the orthonormal basis vectors given by $\frac{1}{\sqrt{n}}e^{2\pi(k)ix}$ evaluated at $n$ regularly spaced points from $[0,1)$. The classic DFT is defined by a basis with a different scaling — namely $\frac{1}{n}$ instead of $\frac{1}{\sqrt{n}}$. This results in the classic DFT having a factor of $\frac{1}{n}$ in the inverse DFT. We choose the convention for the DFT which normalizes in the opposite direction. The basis vectors are just the un-normalized $e^{2\pi(k)ix}$ evaluated at $n$ regularly spaced points from $[0,1)$. We do not want any scaling factors in the relevant inverse DFT because we want to get the coefficients of the Fourier features. of the $n$-dimensional noise vector. By the appropriate212121The reader can verify that the normalization convention we have chosen for the DFT implies $\|\widehat{\bm{\alpha}}\|_{2}^{2}=\frac{1}{n}\|\mathbf{W}_{\mathsf{train}}\|_{2}^{2}$. Parseval’s relation in signal processing, the expected total energy in the feature domain, i.e. $\mathbb{E}\left[\|\widehat{\bm{\alpha}}\|_{2}^{2}\right]$ would be $\sigma^{2}$, and moreover (due to the isotropy of white/independent Gaussian noise), this energy would be equally spread across all $n$ Fourier features in expectation. That is, for every $k\in[n]$ we would have $\mathbb{E}\left[|\widehat{\alpha}_{k}|^{2}\right]=\frac{\sigma^{2}}{n}$.

Now, consider what happens when we include the $\frac{d}{n}-1$ higher-frequency aliases corresponding to each lower frequency component $k\in[n]$. This gives us $d$ Fourier features in total, and we now consider the minimum-$\ell_{2}$-norm interpolator of noise using all $d$ features. The following is what will happen:

In an effort to minimize $\ell_{2}$-norm, the coefficient (absolute) values will be equally divided222222The reader who is familiar with wireless communication theory will be able to relate this to the concept of coherent combining gain. among the $\frac{d}{n}$ aliased features for every realization of the noise samples, i.e. $|\widehat{\alpha}_{k+ln}|=|\widehat{\alpha}_{k}|\text{ for all }l\in\{1,2,\ldots,\frac{d}{n}-1\}$ and for all $k\in[n]$. 2. 2.

For each $k\in[n]$, the expected total contribution from the low frequency feature $k$ and its aliases is now reduced to $\left(\frac{1}{\left(\frac{d}{n}\right)}\right)\cdot\frac{\sigma^{2}}{n}=\frac{\sigma^{2}}{d}$. This results in total $\mathbb{E}\left[\|\widehat{\bm{\alpha}}\|_{2}^{2}\right]=\frac{n}{d}\sigma^{2}$.

For this case, we have zero signal and thus the test MSE for the (whitened) Fourier features is exactly $\|\widehat{\bm{\alpha}}\|_{2}^{2}$. Thus, we have exactly recovered the $\Theta\left(\sigma^{2}\frac{n}{d}\right)$ scaling for the ideal MSE that we bounded in Corollaries 1 and 2. The aliases, when used with minimum-$\ell_{2}$-interpolation of noise, are literally dissipating noise energy, thus directly reducing its potentially harmful effect on generalization in the average232323It is also clear that the average case might be very different than the worst case — a phenomenon intimately connected to the fundamental issue of adversarial examples on neural networks that empirically generalize well. case.

4.1.1 Avoiding signal “bleed"

The problem with $\ell_{2}$-minimizing interpolation is that the above effect of absorbing and dissipating training label energy is generic, whether those labels are signal or noise. Whereas being able to absorb and dissipate training harmful noise energy is a good thing, the same exact effect is harmful when there is true signal. Suppose, as in Equation (14), that the true signal was a constant (thus the only true frequency component is $k=0$). A simple calculation shows that the estimated coefficient of the true function would also be attenuated in exactly the same way, and the absolute value of the coefficient corresponding to the constant feature (i.e. $|\widehat{\alpha}_{0}|=\frac{n}{d}$) decays to [math] as $d\to\infty$. True signal energy, which should ideally be concentrated in the constant feature, is literally bleeding into its aliases in the inference by the $\ell_{2}$-norm-minimizing interpolating solution. This is what we are seeing in the scaling of the test MSE of the $\ell_{2}$-minimizing interpolator for iid Gaussian features (Example 4, shown in Figure 7) as well as the convergence of the test MSE of the minimum-$\ell_{2}$-norm interpolator to the “null risk" proved by Hastie, Tibshirani, Rosset and Montanari [6]. An illustration of this bleeding effect is provided in Figure 4(a), and the realization of this effect on actually recovered signal components for Example 4 is shown in Figure 4(b).

The asymptotic bleeding of signal is a generic issue for whitened features more generally (see [6, Lemma $2$]). This may seem to paint a hopeless picture for the $\ell_{2}$-minimizing interpolator even in the absence of noise – how, then, can it ever work? The key is that we can rescale, and mix, the underlying whitened features to give rise to a transformed feature family, with respect to which we seek an interpolating solution that minimizes the $\ell_{2}$-norm of the coefficients corresponding to these transformed features. The test MSE of such an interpolator will of course be different from the minimum-$\ell_{2}$-norm interpolator using whitened features. The effective difference arises only through the effective rescaling of the whitened features through this transformation: the manifestation of the rescaling can be explicit (the features $\{a(\mathbf{X})_{j}\}_{j=1}^{d}$ can be visibly scaled by weights $w_{j}:=\sqrt{\lambda_{j}}$ for some $\lambda_{j}>0$) or implicit (the eigenvalues of the covariance matrix242424Bartlett, Long, Lugosi and Tsigler [8] present their results through this implicit viewpoint, but their analysis essentially reduces to the explicit viewpoint after a transformation in the underlying geometry. $\bm{\Sigma}$ corresponding to the transformed features $a(\mathbf{X})_{j}$ correspond to the squared weights $\lambda_{j}$).

Consider Example 6 and rescaling $w_{k}>0$ for Fourier feature $f_{k}$. For a set of coefficients $\bm{\alpha}\in\mathbb{R}^{d}$ corresponding to the original whitened features $\{f_{k}\}_{k=0}^{d-1}$, we denote the corresponding coefficients for the rescaled features by $\bm{\beta}\in\mathbb{R}^{d}$. Then, the minimum-$\ell_{2}$-norm interpolator with respect to the rescaled feature family is as below, for any output $\{Y_{j}\}_{j=1}^{n}$:

[TABLE]

This would recover equivalent coefficients $\widehat{\bm{\alpha}}$ for the minimum-weighted-$\ell_{2}$-norm interpolator of the data with weight $\frac{1}{w_{k}}$ corresponding to feature $k$, as below:

[TABLE]

Now, consider the case of a constant signal without noise, i.e. $Y_{j}=1$ for all $j\in[n]$. We already saw that the true signal function, which is $f_{0}(x)$, satisfies $f(x_{j})=1$ for all $j\in[n]$, as does each of its $\left(\frac{d}{n}-1\right)$ aliases $\{f_{ln}(x)\}$ for $l=1,1,\ldots,\frac{d}{n}-1$. Thus, the coefficients of $\widehat{\bm{\alpha}}$ will be the unique linear combination of the aliases, with coefficients represented by $\{\widehat{\alpha}_{ln}\}$, that minimizes the re-weighted $\ell_{2}$-norm subject to the sum of such coefficients being exactly equal to $1$ (to satisfy the interpolation constraint). The special case of whitened features corresponds to $w_{k}=1$ for all $k\in[d]$, and this intuitively results all aliases contributing equally to the recovered signal function. What happens with non-uniform weights: in particular, what happens when $w_{k}$ decreases as a function of frequency $k$? Intuitively, the weighted-$\ell_{2}$-norm objective implies that higher-frequency aliases are penalized more, and thus a solution would favor smaller coefficients $\widehat{\alpha}_{ln}$ for higher integral values of $l$. In fact, Appendix D shows by the principle of matched filtering that the $\ell_{2}$-minimizing coefficients are precisely

[TABLE]

Since the true constant signal is represented by coefficients $\alpha^{*}_{0}=1$ and zero everywhere else, we are particularly interested in the absolute value of $\widehat{\alpha}_{0}$: how much of the true signal component have we preserved? Then, the simple explicit calculation in Appendix D shows that this “survival factor" is essentially252525 This survival factor can also be understood as the outcome of a competition between two functions. The true signal $f_{0}$ that has squared weight $w_{0}^{2}$, and the most attractive orthogonal alias whose squared weight is $\sum_{l=1}^{\frac{d}{n}}w_{ln}^{2}$. The minimum 2-norm interpolator will pick a convex combination of the two by minimizing $\frac{\gamma^{2}}{w_{0}^{2}}+\frac{(1-\gamma)^{2}}{\sum_{l=1}^{\frac{d}{n}}w_{ln}^{2}}$ where $\gamma$ is the survival factor of the true feature. This is minimized by the answer given here for $\gamma=\mathsf{SU}$.

[TABLE]

The inverse of the survival factor $\mathsf{SU}$, after substituting $\lambda_{k}:=w_{k}^{2}$, is very closely related to the first “effective rank" condition introduced by Bartlett, Long, Lugosi and Tsigler to characterize the bias of the minimum-$\ell_{2}$-norm interpolator in [8]. Clearly, the survival factor intimately depends on the relative weights placed on different frequencies, how many frequencies there are in consideration, and how many perfect aliases there are (the number of aliases is inversely proportional to the number of training samples $n$). It is illustrative to rewrite the survival factor $\mathsf{SU}$ as

[TABLE]

Equation (19) is in a form reminiscent of the classic signal-processing “one-pole-filter transfer function”. What matters is the relative weight of the favored feature $w_{0}$ to the combined weight of its competing aliases. As long as it is relatively high, i.e. $w_{0}^{2}\gg\sum_{l=1}^{\frac{d}{n}-1}w_{ln}^{2}$, the true signal will survive. So in particular, if the weights are such that the sum $\sum_{l=1}^{\frac{d}{n}-1}w_{ln}^{2}$ converges even as the number of features grows, the true signal will at least partially survive even as $\frac{d}{n}\to\infty$. Meanwhile, if the sum $\sum_{l=1}^{\frac{d}{n}-1}w_{ln}^{2}$ diverges and does so faster than $w_{0}^{2}$, the signal energy will completely bleed out into the aliases (as happens for the whitened case $w_{k}=1$ for all $k$).

This need for the relative weight on the true features to be high enough relative to their aliases is something that must hold true before any training data has even been collected. In other words, the ability of the 2-norm minimizing interpolator to recover signal is fundamentally restricted. There needs to be a low-dimensional subspace (low frequency signals in our example) that is heavily favored in the weighting, and moreover the true signal needs to be well represented by this subspace. The weights essentially encode an explicit strong prior262626This is in stark contrast to feature selection operators like the Lasso, which select features in a data-dependent manner. that favors low-frequency features.

We can now start to understand the discrepancy between Examples 4 and 5. There is no prior effect favoring in any way the first $500$ features for Example 4. However, by their very nature the features used in Example 5 heavily (implicitly, when the eigenvalue decomposition of $\bm{\Sigma}$ is considered272727In fact, this very case is evaluated in [6, Corollary $1$].) favor the constant feature that best explains the data. This is because the maximal eigenvector of $\Sigma$ is a “virtual feature" that is an average of the $d$ explicit features, i.e. its entries are iid $\mathcal{N}(1,\frac{0.01}{d})$. This better and better approximates the constant feature, the true signal, as $d$ increases – and this improved approximability is the primary explanation for the double descent behavior observed in Figure 3.

In Figure 5, we illustrate how changing the level of the prior weights impacts interpolative solutions using Fourier features for the simple case of a sign function. Here, there is noise in the training data, but the results would look similar even if there were no training noise — the prior weights are primarily fighting the tendency of the interpolator to bleed signal.

4.1.2 Avoiding signal contamination

We have seen that a sufficiently strong prior in a low-dimensional subspace of features avoids the problem of asymptotically bleeding too much of the signal away — as long as the true signal is largely within that subspace. But what happens when some of the true signal is bled away? How does this impact prediction beyond shrinking the true coefficients? Furthermore, the issue of signal bleed does not by itself answer the question of consistency, particularly with the additional presence of noise. How does the strong prior affect fitting of noise – is it still effectively absorbed by the aliases, as we saw when the features were whitened? sTo properly understand this, we need to introduce the idea of “signal contamination.”

Consider Example 6 now with the constant-signal-plus-noise generative model for data:

[TABLE]

The output energy (signal as well as noise) bleeds away from the true signal component corresponding to Fourier feature [math] – but because we are exactly interpolating the output data, the energy has to go somewhere. As a result, all energy that is bled from the true feature will go into the aliased features $\{f_{ln}\}_{l=1}^{d/n-1}$. Each of these features contributes uncorrelated zero-mean unit-variance errors on a test point, scaled by the recovered coefficients $\{\widehat{\alpha}_{ln}\}$. Because they are uncorrelated, their variances add and we can thus define the contamination factor

[TABLE]

Even if there were no noise, the test MSE would be at least $C^{2}$. Consequently, it is important to verify that $C\to 0$ as $(d,n)\to\infty$.

A straightforward calculation (details in Appendix D), again through matched-filtering, reveals that the absolute value of the coefficient on aliased feature $ln$ is directly proportional to the weight $w_{ln}$ and the original true signal strength. Thus contamination (measured as the standard-deviation, rather than the variance in order to have common units), like signal survival, is actually a factor

[TABLE]

for the minimum-weighted-$\ell_{2}$-norm interpolator corresponding to weights $\{w_{k}\}_{k=1}^{d}$. Substituting $\lambda_{k}:=w_{k}^{2}$ results in an error scaling that is very reminiscent of Bartlett, Long, Lugosi and Tsigler’s second effective-rank condition. Thus, we see that the two notions of effective ranks[8] correspond to these factors of survival and contamination, which Bartlett, Long, Lugosi and Tsigler sharply characterize for Gaussian features using random matrix theory. The effective “low-frequency features" there represent directions corresponding to the dominant eigenvalues of the covariance matrix $\bm{\Sigma}$.

There is a tradeoff: while we saw earlier that the weight distribution described here should somewhat favor low-frequency features (dominant eigenvalue directions), it cannot put too little weight on higher-frequency features either. If that happens, the bleeding prevention conditions can be met for a true signal that is in the appropriate low-dimensional subspace. But the noise will give rise to nonvanishing contamination and the variance of the prediction error will not go to [math] as $(n,d)\to\infty$ – then, the minimum-$\ell_{2}$-norm interpolator is inconsistent. Equation (21) tells us that the contamination $C$ will be sufficiently small to ensure consistency iff:

The weights $\{w_{ln}\}_{l\geq 1}$ decay slowly enough so that the sum of squared alias-weights $\sum_{l=1}^{d/n-1}w_{ln}^{2}$ diverges. This means that there is sufficient effective overparameterization to ensure harmless noise fitting. 2. 2.

If the sum of squared alias-weights $\sum_{l=1}^{d/n-1}w_{ln}^{2}$ does not diverge, the term $w_{0}^{2}$ must dominate this sum in the dominator. Then, we also need $w_{0}^{2}\gg\sqrt{\sum_{l=1}^{d/n-1}w_{ln}^{2}}$ so that the denominator dominates the numerator.

Clearly, avoiding non-zero contamination is its own condition, which is not directly implied by avoiding bleeding.

To get consistency, it must be the case that the contamination goes to zero with increasing $n,d$ for everywhere that has true signal as well as an asymptotically complete fraction of the other frequencies. If contamination doesn’t go to zero where the signal is, the test predictions will experience a kind of non-vanishing self-interference from the true signal. If it doesn’t go to zero for most of where the noise is, then that noise in the training samples will still manifest as variance in predictions.

It is instructive to ask whether the above tradeoff in maximizing signal “survival” and minimizing signal “contamination” manifests as a clean bias-variance tradeoff [52]. The issue is that the contamination can arise through signal and/or noise energy. The fraction of contamination that comes from true signal is mathematically a kind of variance that behaves like traditional bias — it is an approximation error that the inference algorithm makes even when there is no noise. The fraction of contamination that comes from noise is indeed a kind of variance that behaves like traditional variance — it would disappear if there were no noise in training data.

4.1.3 A filtering perspective on interpolation

Returning to the case of Fourier features with regularly spaced training points, we can see that given the weightings $w_{i}$ on all the features, we can break the features into cohorts of perfect aliases. All the features are orthogonal (vis-a-vis the test distribution) and because of the regular sampling, each cohort is orthogonal to every other cohort even when restricted to the $n$ sample points. Consequently, we can understand the bleeding within each of the cohorts separately. Moreover, if we assume that the true signal is going to be low-frequency282828This is just for simplicity of exposition and matching the standard machine learning default assumption that all things being equal, we prefer a smoother function to a less smooth function. If the weighting were different, then we could just as well redo this story looking at the highest-weight member of the alias cohort., then we can think about how much the lowest frequency representative of each cohort bleeds. This can be expressed in terms of the survival $0\leq\mathsf{SU}(k^{*})\leq 1$ for that low-frequency feature $k^{*}$ when using the weighted minimum 2-norm interpolator. These $\{\mathsf{SU}(k^{*})\}_{k=0}^{d-1}$ together can be viewed as a filter. This filter tells us how much the act of sampling and estimating attenuates each frequency in the true signal. This attenuation is clearly a kind of “shrinkage.”

With the filtering perspective, we can immediately see that for the minimum-weighted-$\ell_{2}$-norm interpolator to succeed in recovering the signal, the true signal needs to be well approximated by the low-frequency features $\{k^{*}\}$ for which $\mathsf{SU}(k^{*})\approxeq 1$ – otherwise the true pattern will be substantially bled away. We also see that to be able to substantially absorb/dissipate the noise energy (which is going to be spread equally across these $n$ cohorts by the isotropic property of white Gaussian noise), it must be the case that most of the survival coefficients $\{\mathsf{SU}(k^{*})\}_{k=0}^{d-1}$ are quite small — most of the noise energy needs to be bled away. As we tend $(n,d)$ to infinity, we can quantify the required conditions for consistency. In the “continuous time" setting, as $n$ is increasing, the continuous-time frequency (that corresponds to the “fastest" feature) is growing with $n$. So, as long as the maximal value of this frequency $k^{*}$ for which the signal would survive (i.e. $\mathsf{SU}(k^{*})\approxeq 1$) grows sub-linearly292929If this is reminiscent of the conditions discussed when one considers Nadaraya-Watson kernel estimation in nonparametric statistics, this is no coincidence as [53] points out clearly. in $n$, and the set of frequencies for which the signal would “bleed out" (i.e. $\mathbf{SU}(k^{*})\to 0$) is asymptotically $n-o(n)$ frequencies, there is hope of both recovering a low-frequency signal as well as absorbing noise.

On one hand, if $\omega(n)$ of the $\mathsf{SU}$s stay boundedly above [math], then those dimensions of the white noise will clearly not be attenuated as desired, and will show up in our test predictions as a classical kind of prediction variance that is not going to zero. On the other hand, if the true signal is not eventually expressible by low-frequency features whose “survival" coefficients approach $1$, then there is asymptotically non-zero bias in the prediction.

A further nice aspect of the filtering perspective is that it also lets us immediately see that since the relevant Moore-Penrose pseudo-inverse is a linear operator, we can also view it in “time domain.” In machine learning parlance, we could call this the “kernel trick", by which the prediction rule has a direct (and in this case linear) dependence on the labels for the training points. In a traditional signal processing, or wireless communications, perspective, this arises from pulse-shaping filters, or interpolating kernels. A particular set of weights induces both a “survival" filter and an explicit time-domain interpolation function. This is illustrated in Figure 6 for a situation in which we put a substantial prior weight on the low-frequency features. Notice that the low-frequency features survive, and have very little contamination. Meanwhile, the higher-frequecies are attenuated, and though their energy is divided across even higher frequency aliases, the net contamination is also small. The time-domain interpolating kernel looks almost like a classical low-pass-filter, except that it passes through zero at the training point intervals to maintain strict interpolation.

4.1.4 A short comment on Tikhonov regularization and $\ell_{2}$-minimizing interpolation

Hastie, Tibshirani, Rosset and Montanari [6] describe the test MSE of ridge-regularized solutions (with the optimal level of regularization, denoted by $\lambda$) as strictly better than that of the $\ell_{2}$-minimizing interpolator. The discussion above of putting weights on features while we minimize the norm should clearly be evocative of Tikhonov regularization. It turns out that an elementary classical calculation is useful to help us understand what is happening in the new language of bleeding and contamination.

Consider the classical perspective on Tikonov regularization, where here, we will use $\lambda^{2}\Gamma$ as the Tikhonov regularizing matrix to separately call out the ridge-like part $\lambda^{2}>0$ which controls the overall strength of the regularizer and the non-uniformity of the feature weighting that $\Gamma$ represents. As is conventional, let us assume that $\Gamma$ is positive-definite and has a $d\times d$ invertible square-root $\Gamma^{\frac{1}{2}}$ so that $(\Gamma^{\frac{1}{2}})^{\top}\Gamma^{\frac{1}{2}}=\Gamma$. Then, we know that

[TABLE]

In other words, all Tikhonov regularization can be viewed as being a minimum 2-norm interpolating solution for a remixed set of original features combined with the addition of $n$ more special ridge-features that just correspond to a scaled identity — one special feature for each training point. These special “ridge-features” do not predict anything at test-time, but at training time, they do add aliases whose effective expense is controlled by $\frac{1}{\lambda}$. The bigger $\lambda$ is, the cheaper these aliases become, and the more that they bleed signal energy away from other features during minimum 2-norm “interpolative” estimation. Meanwhile, the postmultiplication of $\mathbf{A}_{\mathsf{train}}$ by $\Gamma^{-\frac{1}{2}}$ corresponds to a transformation of the originally given feature family to one that has essentially been premultiplied by $(\Gamma^{-\frac{1}{2}})^{\top}$, causing the new transformed feature family to have covariance $(\Gamma^{-\frac{1}{2}})^{\top}\bm{\Sigma}\Gamma^{-\frac{1}{2}}$. This $\Gamma$ transformation can cause a change in the underlying eigenvalues that is essentially a reweighting.

Consequently, we can understand the Tikhonov-part as essentially being a reweighting of the features. Such a reweighting, by favoring the true parts of the signal and reducing the relative attractiveness of natural aliases, can help control the bleeding if aligned to where the signal actually is. The impact on contamination is indirect and through the same mechanism. Such reweightings can conceivably cause the given feature family’s natural alias structure to be better able to dissipate the noise in unfavored directions. However, the ridge-part is adding additional “aphysical fake features” that are contamination-free by their nature, though they may cause increased bleeding. The contamination-free nature of the ridge-features is coming from the same reason that they cause the prediction to no longer interpolate the training data. Within the context of interpolation, the role of the ridge part is to be a more attractive destination for bled energy than any natural false feature directions while not being attractive at all relative to true feature directions.

This is what tells us immediately that the ridge-part alone cannot prevent attenuation/shrinkage (what we call “bleed”) of signal303030And this is the classical justification for regularizers like Lasso, which stands for “least absolute shrinkage selection operator”!, it can mitigate the additional adverse effect of the attenuated energy manifesting in other features, which we designate as signal contamination. Adding a ridge regularizer cannot help if the main problem is bleeding. This avoidance of contamination, that we have seen could arise either from bled signal or noise, provides a simple explanation for why ridge regularizing is better then $\ell_{2}$-minimizing interpolation [6]. It is also a nice high-level explanation for the double-descent behavior of the Tikhonov regularizer in [46].

5 Interpolation in the noisy sparse linear model

From the results obtained so far, we can think of any successful interpolator as recovering the signal in the noiseless high-dimensional regime, and then fitting noise harmlessly as in Corollary 2. We just saw that the $\ell_{2}$-interpolator is generically extremely poor at identifying the correct signal $\bm{\alpha}^{*}$ even in the noiseless setting – this was fundamentally due to the issue of the signal “bleeding out" across too many features, leading to practically a zero-fit as the overparameterization went to infinity. If the signal $\bm{\alpha}^{*}$ is dense, we know that there is information-theoretically no hope of identifying it [41, 42]. On the other hand, if there is underlying sparsity in the signal we know that recovery is possible with estimators that are fundamentally sparsity-seeking, i.e. estimators obtained through procedures that induce parsimony in the coefficients of $\widehat{\bm{\alpha}}$. Traditionally, the estimators used do not interpolate in the noisy regime – they regularize to essentially denoise the output signal313131(Common examples are the Lasso and orthogonal matching pursuit (OMP) with an early stopping rule.) . In other words, fitting noise has always been thought of as harmful. We now ask whether this is actually the case by analyzing estimators that interpolate in the presence of noise. First, we define the sparse linear model.

Definition 3** (Noisy sparse linear model).**

For any $k\geq 0,\sigma>0$, the $(k,\sigma)$-sparse whitened linear model describes output that is generated as

[TABLE]

where $W\sim\mathcal{N}(0,\sigma^{2})$ and $\bm{\alpha}^{*}$ is $k$-sparse, i.e. $\|\bm{\alpha}^{*}\|_{0}\leq k$. We denote the support of $\bm{\alpha}^{*}$ by $S^{*}:=\mathsf{supp}(\bm{\alpha}^{*})$ and note that $|S^{*}|\leq k$. For ease of exposition, we also assume $\mathbb{E}\left[\mathbf{a}(\mathbf{X})\mathbf{a}(\mathbf{X})^{\top}\right]=\mathbf{I}_{d}$, i.e. the input features are whitened.

5.1 Sparsity-seeking methods run to completion

A starting choice for a reasonable practical interpolator in the sparse regime might be an estimator meant for the noiseless sparse linear model to fit the signal as well as noise. Two examples of such estimators are below:

Orthogonal matching pursuit (OMP) to completion. 2. 2.

Basis pursuit (BP): ${\arg\min}\|\bm{\alpha}\|_{1}\text{ subject to Equation~{}\eqref{eq:interpolatingsoln}.}$

The above methods are extensively analyzed for noiseless sparse linear models [45], i.e. the $(k,0)$-sparse linear model as defined in Definition 3. While they are also known to successfully recover signal in the high-signal/zero-noise limit, less is known about their performance in the presence of substantial noise. The test MSE of these interpolators as a function of $d$ in this case is shown in Figure 7. We observe that OMP run to completion and BP are significantly better than the minimum $\ell_{2}$ norm interpolator. In fact, they appear only an extra noise multiple of test MSE over the ideal interpolator. The vast improvement of these estimators over the minimum $\ell_{2}$-norm interpolator is not surprising on some level, as at the very least these estimators will preserve signal. The estimator that uses BP in its first step, and then thresholds the coefficients that are below (half the) minimum absolute value of the non-zero coefficients of $\bm{\alpha}^{*}$, has been shown to recover the sign pattern, or the true support, of $\bm{\alpha}^{*}$ exactly [54]. Indeed, this implies that after the first step, the true signal components (corresponding to entries in $\mathsf{supp}(\bm{\alpha}^{*})$) are guaranteed to be preserved by BP even in the presence of low-enough levels of noise. This means that the signal is preserved, or “survives" to use our language, with these sparsity-seeking interpolators in the presence of noise. Similarly, one can show that the first $k$ steps of OMP run to completion will indeed select the correct features corresponding to $\mathsf{supp}(\bm{\alpha}^{*})$ with high probability; for more details about OMP, see Appendix A.3.

What is not immediately clear about BP and OMP run to completion is their effect on fitting the noise itself – does it overfit terribly, or harmlessly (like in Corollary 2)? This is directly connected to the “signal contamination" factor discussed in Section 4, and is not directly answered by existing analysis. For example, [54] only guarantees that spurious coefficients are at most half the minimum non-zero entry of $\bm{\alpha}^{*}$, which is generally a constant. To get a clearer picture of what may happen to purely sparsity-seeking interpolators, we isolate the effect of sparsity-seeking interpolation on noise. Consider the special case of zero signal, i.e. let $\bm{\alpha}^{*}=\mathbf{0}$. In this case, for whitened feature families the test MSE of an estimator $\widehat{\bm{\alpha}}$ is simply $\|\widehat{\bm{\alpha}}\|_{2}^{2}$, and the estimator $\widehat{\bm{\alpha}}$ is in fact fitting pure noise. Figure 8 shows the scaling of the test MSE in the zero-signal case as a function of $n$ and $d$ for OMP and BP. The ideal test MSE, i.e., the fundamental price of interpolation of noise, is also plotted for reference. We make these plots for three high-dimensional regimes:

$n$ fixed, and $d\geq n$ growing. In this regime, we ideally want the test MSE to decay to [math] as $d\to\infty$. This regime is plotted in Figure 8(a) for $n=1000$, and 8(b) for $n=50$. 2. 2.

$d=n^{2}$ and $n$ growing. We ideally want consistency in the sense that we want the test MSE to decay to [math] as $(n,d)\to\infty$. This regime is plotted in Figure 8(c). 3. 3.

$d=e^{n}$ and $n$ growing. As before, we want the test MSE to to decay to [math] as $(n,d)\to\infty$. This regime is plotted in Figure 8(d).

In all the regimes, Figure 8 shows us that the test MSE of the sparsity-seeking interpolators does slightly decrease with $(n,d)$, but extremely slowly and negligibly in comparison to the ideal test MSE. The issue is that these interpolators are fundamentally parsimonious: they use $\mathcal{O}(n)$ features to fit the noise. While this property was desirable for signal recovery, it is not that desirable for fitting noise as harmlessly as possible. We define such an “overly parsimonious interpolating operator" broadly below.

Definition 4**.**

For a fixed training data matrix $\mathbf{A}_{\mathsf{train}}$, consider any interpolating operator $\widehat{\bm{\alpha}}:=\widehat{\alpha}(\mathbf{Y})$. Let $S_{\mathsf{top},n}:=\{s_{j}\}_{j=1}^{n}\subset[d]$ represent the indices for the top $n$ absolute values of coefficients of $\widehat{\bm{\alpha}}$. Then, define truncated vector $\widehat{\bm{\alpha}}_{\mathsf{trunc}}$ such that

[TABLE]

Then, for some constant $\beta\in(0,1]$, a $\beta$-parsimonious interpolating operator satisfies

[TABLE]

for all vectors $\mathbf{Y}\in\mathbb{R}^{n}$. For the special case $\beta=1$, we will call this a hard-sparse interpolating operator.

Definition 4 applies to interpolating operators that are a function of the output $\mathbf{Y}_{\mathsf{train}}$. For a given data matrix, a parsimonious interpolating operator selects exactly $n$ features that capture a constant fraction of the energy of any output. (The locations of the $n$ features that are selected will, of course, depend on the output.) For the special case of $\beta=1$, we will see that hard-sparse interpolators include OMP and the minimum-$\ell_{1}$-norm interpolator (BP). Because a non-vanishing fraction (represented by the parameter $\beta>0$) of the noise energy is concentrated in only $n$ features, the effective overparameterization benefit of ideal interpolation of noise (that we saw in Corollary 2) can no longer be realized.

We state the main result of this section for a random whitened feature matrix $\mathbf{A}_{\mathsf{train}}$ satisfying one out of sub-Gaussianity (Assumption 2) or Gaussianity (Assumption 3).

Theorem 2**.**

Consider any interpolating solution $\widehat{\bm{\alpha}}$ obtained by a $\beta$-parsimonious interpolating operator (as defined in Definition 4) with constant $\beta\in(0,1]$. Then, when applied to any random whitened feature matrix $\mathbf{A}_{\mathsf{train}}$ satisfying:

Gaussianity (Assumption 3), there exists an instance of the $(k,\sigma^{2})$-sparse linear model for any $k\geq 0$ for which the test MSE

[TABLE]

for any $\delta>0$, with probability at least $(1-e^{-n\ln\left(\frac{d}{n}\right)}-e^{-n\delta^{2}/2})$ over realizations of feature matrix $\mathbf{A}_{\mathsf{train}}$ and noise $\mathbf{W}_{\mathsf{train}}$. 2. 2.

sub-Gaussianity of rows (Assumption 2) with parameter $K>0$, there exists an instance of the $(k,\sigma^{2})$-sparse linear model for which the test MSE

[TABLE]

for any $\delta>0$, with probability at least $(1-e^{-n\ln\left(\frac{d}{n}\right)}-e^{-n\delta^{2}/2})$ over realizations of feature matrix $\mathbf{A}_{\mathsf{train}}$ and noise $\mathbf{W}_{\mathsf{train}}$. Here, constant $C^{\prime\prime}_{K}>0$ depends only on the upper bound on the sub-Gaussian parameter, $K$.

Theorem 2 should be thought of as a negative result for the applicability of parsimonious interpolators meant for the noiseless setting in additionally fitting noise, even when the setting is heavily overparameterized – that is, we are no longer enjoying harmless interpolation of noise. Consider the following scalings of $d$ with respect to $n$:

$d=\gamma n$ for some constant $\gamma>1$. In this case, Equations (27) and (28) give us $\|\widehat{\bm{\alpha}}\|_{2}^{2}=\omega(\sigma^{2})$, and the test MSE does not go to [math] as $n\to\infty$. 2. 2.

$d=n^{q}$ for some $q>1$. In this case, Equations (27) and (28) give us $\|\widehat{\bm{\alpha}}\|_{2}^{2}=\omega\left(\frac{\sigma^{2}}{(q-1)\ln n}\right)$ which goes to [math] as $n\to\infty$, but at an extremely slow logarithmic rate. 3. 3.

$d=e^{\gamma n}$ for some $\gamma>0$. This is an extremely overparameterized regime. In this case, Equations (27) and (28) give us $\|\widehat{\bm{\alpha}}\|_{2}^{2}=\omega\left(\frac{\sigma^{2}}{\gamma n-\ln n}\right)$ which is a much faster rate. However, in this exponentially overparameterized regime it is well known that successful signal recovery is impossible even in the absence of noise.

Putting these conclusions together, Theorem 2 suggests that while consistency of parsimonious interpolators might be possible in polynomially high-dimensional regimes – it would be at an extremely slow logarithmic rate.

Theorem 2 holds for a broad class of sparsity-seeking interpolators that successfully recover signal in the absence of noise (in the polynomially high-dimensional regime). The following results hold for OMP run to completion, and BP – both of which satisfy $1$-parsimonious interpolation on any output.

Corollary 3**.**

When the random whitened feature matrix $\mathbf{A}_{\mathsf{train}}$ satisfies one out of Gaussianity (Assumptions 3) or sub-Gaussianity (Assumption 2), the interpolator formed by OMP run to completion incurs test MSE

[TABLE]

for some constant $C>0$ with high probability for at least one instance of the $(k,\sigma)$-sparse linear model.

Corollary 3 is trivial because, by nature, OMP run to completion selects exactly $n$ features and stops (see Appendix A.3 for details.) Thus, the OMP interpolator is always $n$-hard-sparse, thus $1$-parsimonious according to Definition 4.

Corollary 4**.**

When the random whitened feature matrix $\mathbf{A}_{\mathsf{train}}$ satisfies one out of Gaussianity (Assumptions 3) or sub-Gaussianity (Assumption 2), the minimum-$\ell_{1}$-norm interpolator (BP) incurs test MSE

[TABLE]

for some constant $C>0$ with high probability for at least one instance of the $(k,\sigma)$-sparse linear model.

Corollary 4 is proved in Appendix A.2. The result follows trivially from Theorem 2 once we know that the minimum-$\ell_{1}$-norm interpolator also selects $n$ features. While this is less immediately obvious than OMP, it is a classical result323232Also referenced in Chen, Donoho and Saunder’s classic survey [45]. The name “basis pursuit” is actually a give-away for this property, as we know that a basis in $\mathbb{R}^{n}$ cannot contain more than $n$ elements., that follows from the linear program (LP) formulation of BP, that the $\ell_{1}$-minimizing interpolator selects at most $n$ features, and will generically select exactly $n$ features. We reproduce this argument in full in Appendix A.2 for completeness.

We close this section with the proof of Theorem 2.

Proof.

For any $k>0$, we consider the instance of the $(k,\sigma)$-sparse linear model for which there is zero signal, i.e. $\bm{\alpha}^{*}=\mathbf{0}$. In this case, the output is pure noise, i.e. $\mathbf{Y}_{\mathsf{train}}=\mathbf{W}_{\mathsf{train}}$ and the interpolator operates on pure noise as $\widehat{\bm{\alpha}}:=\alpha(\mathbf{Y}_{\mathsf{train}})=\alpha(\mathbf{W}_{\mathsf{train}})$. Further, recall that the test MSE on whitened features for any estimator $\widehat{\bm{\alpha}}$ is defined as

[TABLE]

Recall that $\widehat{\bm{\alpha}}_{\mathsf{trunc}}$ is the truncated version of the interpolator $\widehat{\bm{\alpha}}$ as defined in Equation (25). Let $S=\mathsf{supp}(\widehat{\bm{\alpha}}_{\mathsf{trunc}})$ denote the support of the truncated vector $\widehat{\bm{\alpha}}_{\mathsf{trunc}}$. Recall, that by definition, $|S|=n$. More generally, the actual composition of the $n$ elements in $S$ will depend both on the realizations of the random matrix $\mathbf{A}_{\mathsf{train}}$ and the noise $\mathbf{W}_{\mathsf{train}}$.

Let $\mathbf{A}_{\mathsf{train}}(S)$ be the $n\times n$ matrix with columns sub-sampled from the set $S$ of size $n$, and denote $\mathbf{W}_{\mathsf{train}}^{\prime}:=\mathbf{A}_{\mathsf{train}}\widehat{\bm{\alpha}}_{\mathsf{trunc}}=\mathbf{A}_{\mathsf{train}}(S)\widehat{\bm{\alpha}}_{\mathsf{trunc}}$. Assuming that $\mathbf{A}_{\mathsf{train}}(S)$ is invertible (which is always almost surely true for a random matrix), we note that

[TABLE]

and so we have

[TABLE]

Thus, we need to prove a lower bound on the maximal eigenvalue $\lambda_{max}\left(\mathbf{A}_{\mathsf{train}}(S^{\prime})^{\top}\mathbf{A}_{\mathsf{train}}(S^{\prime})\right)$ that will hold, point-wise, for all $S^{\prime}\subset[d],|S^{\prime}|=n$. We state and prove the following intermediate lemma.

Lemma 2**.**

For matrix $\mathbf{A}_{\mathsf{train}}$ satisfying:

Assumption 3 (Gaussian features), we have

[TABLE]

with probability greater than or equal to $(1-e^{-n\ln\left(\frac{d}{n}\right)})$. 2. 2.

Assumption 2 (Sub-Gaussianity), we have

[TABLE]

with probability greater than or equal to $(1-e^{-n\ln\left(\frac{d}{n}\right)})$, where parameter $C_{K},c_{K}>0$ are positive constants that depend on the sub-Gaussian parameter $K$.

Notice that Lemma 2 directly implies the statement of Theorem 2. This is because with probability greater than or equal to $(1-e^{-n\ln\left(\frac{d}{n}\right)})$, we then have for any subset $|S^{\prime}|=n$,

[TABLE]

for positive constants $C,C^{\prime},C^{\prime\prime}$ (we have dropped the subscript dependence on $K$ for the sub-Gaussian case for convenience). Now, recall that $\mathbf{W}_{\mathsf{train}}^{\prime}=\mathbf{A}_{\mathsf{train}}\widehat{\bm{\alpha}}_{\mathsf{trunc}}$, and from Equation (26) in Definition 4, we get

[TABLE]

Also recall from the lower tail bound on chi-squared random variables that $\|\mathbf{W}_{\mathsf{train}}\|_{2}^{2}\geq n\sigma^{2}(1-\delta)$ with probability greater than or equal to $(1-e^{-\frac{n\delta^{2}}{8}})$. Putting these facts together, we get

[TABLE]

Finally, we note that truncation to [math] only decreases the $\ell_{2}$-norm, and thus $\|\widehat{\bm{\alpha}}\|_{2}^{2}\geq\|\widehat{\bm{\alpha}}_{\mathsf{trunc}}\|_{2}^{2}$. This gives us

[TABLE]

which is precisely the statement in Theorem 2. Note that the overall probability of this statement is greater than or equal to $\left(1-e^{-\frac{n\delta^{2}}{8}}-e^{-n\ln\left(\frac{d}{n}\right)}\right)$.

It only remains to prove Lemma 2, which we do below.

Proof.

Under Gaussianity (Assumption 3), we use Lemma 1 which we introduced in the proof of Corollary 1. For every $S^{\prime}\subset[d],|S^{\prime}|=n$, a direct substitution of of Lemma 1 applied to the maximum singular value of matrix $\mathbf{A}_{\mathsf{train}}(S^{\prime})$ with $t:=t_{0,1}=\sqrt{4n\left(\ln\left(\frac{d}{n}\right)+1\right)}$ gives us

[TABLE]

Similarly, under Assumption 2, we use Lemma 3 stated in Appendix A.1. For every $S^{\prime}\subset[d],|S^{\prime}|=n$, a direct substitution of of Lemma 3 applied to the maximum singular value of matrix $\mathbf{A}_{\mathsf{train}}(S^{\prime})$ with $t=t_{0,2}:=\sqrt{\frac{2n\left(\ln\left(\frac{d}{n}\right)+1\right)}{c_{K}}}$ gives us

[TABLE]

For convenience, we unify the rest of the argument for both assumptions, denoting $C:=\{1,C_{K}\}$, $c:=\{\frac{1}{2},c_{K}\}$ and $t_{0}:=\{t_{0,1},t_{0,2}\}$ for the Gaussian case and sub-Gaussian case respectively. Applying the union bound for all $S^{\prime}\subset[d],|S^{\prime}|=n$ gives us

[TABLE]

where we have substituted $t_{0}:=\sqrt{\frac{2n\left(\ln\left(\frac{d}{n}\right)+1\right)}{c}}$ and inequality $(\mathsf{i})$ follows from Fact 1 in Appendix C. This completes the proof for both the Gaussian case and the sub-Gaussian case. ∎

Notice that the union-bound used here is quite likely a bit loose. But we know that the maximum of independent random variables does not behave far from what the union-bound predicts, and so tightening here would require exploiting non-independence. ∎

5.2 Order-optimal hybrid methods

Can we construct an optimal interpolation scheme in the $(k,\sigma)$-noisy sparse linear model for any $(k,\sigma)$? We define optimality in the sense of order-optimality, i.e. the optimal scaling for test MSE as a function of $(d,n)$ that holds for all instances333333Statisticans commonly call this minimax-optimality. in the $(k,\sigma)$-noisy sparse linear model. Clearly, any order-optimal interpolator needs to satisfy the following two conditions:

The test MSE due to the presence of signal should be proportional to $\sigma^{2}\cdot\frac{k\ln\left(\frac{d}{k}\right)}{n}$. This is the best possible scaling in test MSE that any estimator can achieve in the noisy sparse linear model, whether or not such an estimator interpolates. 2. 2.

The test MSE due to the presence of noise, even when the signal is absent, should be proportional to $\sigma^{2}\cdot\frac{n}{d}$. As we saw in Corollaries 1 and 2, this is the fundamental price of interpolation of noise.

Let us consider the requirements of an interpolator to meet the above two conditions for order-optimality. For one, it should be able to uniquely identify the signal in the absence of noise; two, it needs to fit noise harmlessly across $\omega(d)$ features. We saw in Section 4 that the minimum-$\ell_{2}$-norm interpolator fails at signal identifiability because it tends to bleed signal out across multiple features, rather than concentrating it in a few entries. On the other hand, in Section 5.1 we saw that purely sparsity-seeking interpolators are poor choices to fit noise harmlessly because of their parsimony – they only spread noise among $n$ features. Thus, the goals for signal identifiability and harmless noise-fitting are fundamentally at odds with one another in the overparameterized regime with underlying sparsity.

This tradeoff suggests that the “optimal" interpolator should use different procedures for fitting the part of the output that is signal, and the part of the output that is noise. It suggests the application of a hybrid scheme, described in the $k$-sparse regime below. Recall that we have assumed whitened feature families, i.e. $\bm{\Sigma}=\mathbf{I}_{d}$, for ease of exposition343434This is primarily to be able to state our corollaries for sparse signal recovery out-of-the-box. However, more general versions of these results will hold for unwhitened feature families as well. .

Definition 5**.**

Corresponding to any estimator $\widehat{\bm{\alpha}}_{1}$ of $\bm{\alpha}^{*}$ for the $(k,\sigma$)-whitened sparse linear model (note that this estimator need not interpolate), we can define a two-step hybrid interpolator as below:

Compute the residual $\mathbf{W}_{\mathsf{train}}^{\prime}=\mathbf{Y}_{\mathsf{train}}-\mathbf{A}_{\mathsf{train}}\widehat{\bm{\alpha}}_{1}$. 2. 2.

Return \widehat{\bm{\alpha}}_{\mathsf{hybrid}}:={\arg\min}\|\bm{\alpha}-\widehat{\bm{\alpha}}_{1}\|_{2}\text{ subject to \bm{\alpha} satisfying Equation~{}\eqref{eq:interpolatingsoln} }. Clearly, an equivalent characterization is $\widehat{\bm{\alpha}}_{\mathsf{hybrid}}=\widehat{\bm{\alpha}}_{1}+\widehat{\bm{\Delta}}$, where $\widehat{\bm{\Delta}}:={\arg\min}\|\bm{\Delta}\|_{2}\text{ subject to }\mathbf{A}_{\mathsf{train}}\bm{\Delta}=\mathbf{W}_{\mathsf{train}}^{\prime}$.

Observe that the feasibility constraint is just a rewriting of Equation (1), ensuring that the estimator $\widehat{\bm{\alpha}}_{\mathsf{hybrid}}$ interpolates the data. The guarantee that such a hybrid scheme achieves on test MSE is stated in the following proposition.

Proposition 1**.**

Denote the estimation error and prediction error of the estimator $\widehat{\bm{\alpha}}_{1}$ by

[TABLE]

Then, for any estimator $\widehat{\bm{\alpha}}_{1}$, the hybrid estimator $\widehat{\bm{\alpha}}_{\mathsf{hybrid}}$ has test MSE

[TABLE]

Proposition 1 shows a natural split in the error of such hybrid interpolators in terms of two quantities: the error that arises from signal recovery, and the error that arises from fitting noise. The proof of Proposition 1 follows simply from the insights already presented and is given in Appendix A.3 for completeness.

Using the ideas from the proof of Corollary 2, one can then show that the performance of a hybrid interpolator fundamentally depends on the ability of the sparse signal estimator to estimate the signal, plus the excess error incurred by the ideal interpolator of pure noise, i.e. the ideal test MSE arising from fitting noise. For the rest of this section, we focus on the special case where the entries of $\mathbf{A}_{\mathsf{train}}$ are drawn from the standard normal distribution, i.e. $\mathbf{A}_{\mathsf{train}}$ satisfies Assumption 3 together with whiteness. The literature on sparse signal recovery is rich, and we can interpret the guarantee provided by Proposition 1 for a variety of choices of the first-step estimator $\widehat{\bm{\alpha}}_{1}$. A summary of these results is contained in Table 1.

First, we consider estimators $\widehat{\bm{\alpha}}_{1}$ that are optimal in their scaling with respect to $(k,\sigma^{2},n,d)$ in the sparse regime. From here on, we will call these order-optimal estimators.

Corollary 5**.**

Consider the feature matrix $\mathbf{A}_{\mathsf{train}}$ with iid standard Gaussian entries, and any estimator $\widehat{\bm{\alpha}}_{1}$ that is order-optimal in both estimation error and prediction error, i.e. there exist universal constants $C,C^{\prime}>0$ such that

[TABLE]

with high probability (over the randomness of both the noise $\mathbf{W}_{\mathsf{train}}$ and the randomness in the whitened feature matrix $\mathbf{A}_{\mathsf{train}}$). Then, provided that $d>4n$, the hybrid estimator $\widehat{\bm{\alpha}}_{\mathsf{hybrid}}$ that is based on estimator $\widehat{\bm{\alpha}}_{1}$ gives us test MSE

[TABLE]

for universal constants $C^{\prime\prime},C^{\prime\prime\prime}>0$. Such a hybrid interpolator is order-optimal among all interpolating solutions.

To see that the rate in Equation (31) is order-optimal, note that any interpolating solution would need to incur error at least $\sigma^{2}\frac{n}{(\sqrt{d}+2\sqrt{n})^{2}}$ due to the effect of fitting pure noise (from the lower bound part of Corollary 2). On the other hand, we know that any sparse-signal estimator, regardless of whether it interpolates or not, would need to incur error at least $\sigma^{2}\frac{k\ln\left(\frac{d}{k}\right)}{n}$. Thus, the test MSE of any interpolating solution on at least one instance in the sparse regime would have to be at least

[TABLE]

which exactly matches the rate in Equation (31) upto constants.

One example of such an order-optimal estimator is the SLOPE estimator which was recently analyzed [55]. If one is willing to tolerate a slightly slower rate of $\mathcal{O}\left(\sigma^{2}\cdot\frac{k\ln d}{n}\right)$ for signal recovery, several other estimators can be used. We state the recovery guarantee (informally) for two choices below. Formal statements are proved in Appendix A.3.

Corollary 6**.**

Consider the estimators $\widehat{\bm{\alpha}}_{1}:=\widehat{\bm{\alpha}}_{\mathsf{Las.},1}$ and $\widehat{\bm{\alpha}}_{1}:=\widehat{\bm{\alpha}}_{\mathsf{OMP},1}$ obtained by Lasso and OMP respectively for suitable choices of regularizer and stopping rule353535Details in Appendix A.3. respectively. Assuming lower bounds on number of samples $n$ and signal strength respectively, the hybrid estimator based on either of these estimators then gives us test MSE

[TABLE]

for universal constants $C^{\prime\prime},C^{\prime\prime\prime}>0$.

Corollaries 5 and 6 are proved in Appendix A.3 with explicit analysis of Lasso and OMP for completeness. The signal recovery guarantees are slightly different for the Lagrangian Lasso and OMP - for the former, we restate a direct bound [56] on the $\ell_{2}$-norm of the signal recovery error, which scales as given in the informal equation above. For the latter, our error bound derives from a much stronger variable selection guarantee [57] which implicitly requires the SNR to scale as $\omega(\ln d)$. For both, it is worth paying attention to the conditions strictly required for successful recovery – for the Lagrangian Lasso, the condition manifests as a lower bound on the number of samples $n$ as a function of $(k,d)$. For OMP, the condition manifests as a lower bound on the absolute values of non-zero coefficients $\{\alpha^{*}_{j}\},j\in\mathsf{supp}(\bm{\alpha}^{*})$. For further details on these conditions, see Appendix A.3.

The examples of estimators provided so far for the first step of the hybrid interpolator require knowledge of the noise variance $\sigma^{2}$ – they use it either to regularize appropriately, or to define an appropriate stopping rule. While one could always estimate this parameter through cross-validation, there do exist successful estimators that do not have access to this side information. One example is provided below.

Corollary 7**.**

Consider the estimator $\widehat{\bm{\alpha}}_{1}:=\widehat{\bm{\alpha}}_{\mathsf{Sq.rt.Las.},1}$ which is based on the square-root-Lasso for suitable choice of regularizer that does not depend on the noise variance $\sigma^{2}$ (details in Appendix A.3). The hybrid estimator then gives us test MSE

[TABLE]

for universal constants $C^{\prime\prime},C^{\prime\prime\prime}>0$.

Regardless of which sparse signal estimator is used in the first step of the hybrid interpolator (Definition 5), the qualitative story is the same: there is a tradeoff in how to overparameterize, i.e. how to set $d>>n$. As we increase the number of features in the family, the error arising from fitting noise goes down as $1/d$ - but these “fake features" also have the potential to be falsely discovered363636This same idea inspired recent work on explicitly constructing fake “knockoff” features to draw away energy from features that have the potential to be falsely discovered [58]., thus driving up the signal recovery error rate logarithmically in $d$. This logarithmic-linear tradeoff still ensures that the best test error is achieved when we sizeably overparameterize, even if not at $d=\infty$ like if we were only fitting noise. Figure 7 considers an example with iid Gaussian design and true sparsity level $500$ for various ranges of noise variance: $\sigma^{2}=10^{-4},10^{-2}$ and $10^{-1}$. For all these ranges, we observe that hybrid interpolators’ test MSE closely track the ideal MSE.

6 Conclusions for high-dimensional generative model

Our results, together with other work in the last year [6, 7, 8, 9], provide a significant understanding of the ramifications of selecting interpolating solutions of noisy data generated from a high-dimensional, or overparameterized linear model, i.e. where the number of features used far exceeds the number of samples of training data. Key takeaways are summarized below:

While “denoising" the output is always strictly preferred to simply interpolating both signal and noise, the additional price of this interpolation on noise can be minimal when the regime is extremely overparameterized. The price decays to [math] as the level of overparameterization goes to infinity. We can now intuitively see that this is a consequence of the ability of aliased features to absorb and dissipate training noise energy — figuratively spreading it out over a much larger bandwidth. 2. 2.

Realizing this potential for “harmless interpolation of noise” critically depends on the number of such aliased features, as shown for Gaussian feature families in [8] for the minimum $\ell_{2}$-norm interpolator, and this in fact holds completely generally even for ideal interpolation. Overparameterization is necessary, for harmless interpolation. 3. 3.

The $\ell_{2}$-minimizing interpolator of noisy data with whitened features successfully absorbs noise energy for maximal harmlessness, but also dissipates signal across multiple aliases, as also shown in [6]. To prevent this “bleeding of signal", there needs to be a strong implicit bias towards a small, a-priori known set of pre-determined feature directions. Moreover, as mentioned in [8], to be successful, we should only care about measuring generalization error along these pre-determined feature directions — i.e. the true signal needs to be in this known subspace to be reliably learned by the $\ell_{2}$-min interpolator. 4. 4.

The $\ell_{1}$-minimizing interpolator of noisy data suffers from the reverse problem: it more successfully preserves signal (when indeed the true signal is sparsely representable), but does not absorb training noise harmlessly. 5. 5.

Interpolators in the sparse linear model satisfying the diametrically opposite goals of sparse signal preservation and noise energy absorption exist, and can be understood as having a two-step, “hybrid" nature. They first fit the signal as best they can and then interpolate everything in a way that spreads noise out. This gives asymptotic rates that are, in an order sense, the best that could be hoped for.

These takeaways paint a more measured picture for the practical use of interpolation in noisy high-dimensional linear regression. While it is intellectually interesting that interpolation need not be harmful, it is practically always suboptimal to regularization with denoising. This was pointed out in [6, Sections $6$ and $7$] for the $\ell_{2}$-minimizing interpolator vis-a-vis Tikhonov regularization at the optimal level (and this optimal level can be estimated using cross-validation). Our discussion of Tikhonov and ridge regularization shows that this is because ridge regression can be interpreted as a further level of overparameterization, except using features that manifest as [math] during test-time. The best possible interpolating estimators are an additional $\mathcal{O}\left(\sigma^{2}\frac{n}{d}\right)$ worse in test MSE than their optimal denoising counterparts.

Moreover, because of the fundamental tradeoff between signal preservation and noise absorption, constructed consistent estimators in a two-step procedure by building on optimal regularizers to further fit noise – a second step that indeed appears rather artificial and unnecessary. A conceivable benefit to interpolating solutions was their lack of dependence on prior knowledge of the noise variance $\sigma^{2}$; however, regularizing estimators are also now known to work in the absence of this knowledge – either by modifying the optimization objective to make it SNR-invariant [59], or by first estimating the SNR [60].

6.1 Future Directions

The benefits of overparameterization, however, deserve considerable future attention and scrutiny – regardless of whether the solutions interpolate or regularize. The missing ingredient in our discussion above is the potential for model misspecification, i.e. data that can only be approximated by the high-dimensional linear model. (The example 5 in Section 4 is just a caricature in that direction.) In the empirical double descent papers [37, 38], the second descent of the test MSE in the overparameterized regime is strictly better than the first descent in the underparameterized regime. From our concrete understanding of the correctly specified high-dimensional regime, it is clear that this is only possible in the presence of an approximation-theoretic benefit of adding more features into the model. Analyzing this benefit is very challenging generally, but in the special case when the features are independent, [6, Theorem $5$] provides asymptotic minimal scalings for the approximation-theoretic benefit under which a double descent curve can be obtained with the $\ell_{2}$-minimizing interpolator. Even here, the global minimum of the empirical risk is not at infinite overparameterization. Very recetly, Mei and Montanari [46] provided a promising result: they theoretically recovered the second descent of the double-descent curve for the Tikhonov-regularized solution (note that this is not an interpolator) on random Fourier features, which are clearly not independent. More work remains to be done to characterize necessary and sufficient conditions for a double descent in the empirical risk – such conditions include the generative model for the data, the choice of featurization, the choice of estimator, and whether this estimator should even interpolate.

Finally, as also pointed out by [8], we are far from understanding the ramifications of overparameterization on generalization for the original case study that motivated recent interest in the overparameterized regime: deep neural networks with many more parameters than training points [11, 4]. In most practically used neural networks, the overparameterization includes bothdepth rather than width, resulting in the presence of substantial non-linearity – how this affects generalization, even in this intuitive and simple picture of signal preservation and noise absorption – remains unclear. As a starting point to investigate this case, [6, Section $8$] provides a model that interpolates between full linearity and non-linearity.

Quantifiable ramifications of interpolation and overparameterization on non-quadratic loss functions (such as those that appear in classification problems) would also be very interesting to understand, but the investigation here suggests that similar effects (bleeding, noise absorption, etc.) should exist there as well since the core intuition underlying them is that they are consequences of aliasing. The signal-processing perspective here actually suggests a pair of families of meta-conjectures that can help drive forward our understanding. On the one hand, we can use the lens of regularly-sampled training data and Fourier features to understand phenomena of interest whenever we seek to understand the significantly overparameterized regime in the context of inference algorithms that are minimum 2-norm in nature. The resulting conjectures can then be validated in Gaussian models and beyond. But perhaps even more interesting is the reverse direction — our lens suggests that the many phenomena and techniques developed over the decades in the regularly-sampled Fourier world should have counterparts in the world of overparameterized learning.

Acknowledgments

We thank Peter L Bartlett for insightful initial discussions, and Aaditya Ramdas for thoughtful questions about an earlier version of this work. We would like to acknowledge the staff of EECS16A and EECS16B at Berkeley for in part inspiring the initial ideas, as well as the staff of EECS189.

We thank the anonymous reviewers of IEEE ISIT 2019 for useful feedback, and both the information theory community and the community at the Simons Institute Foundations of Deep Learning program, Summer 2019, for several stimulating discussions that improved the presentation of this paper – especially Andrew Barron, Mikhail Belkin, Shai Ben-David, Meir Feder, Suriya Gunasekar, Daniel Hsu, Tara Javidi, Ashwin Pananjady, Shlomo Shamai, Nathan Srebro, and Ram Zamir.

Last but not least, we acknowledge the support of the ML4Wireless center member companies and NSF grants AST-144078 and ECCS-1343398.

Appendix A Supplemental proofs

We start by formally stating the matrix concentration results that are used in the proof of Corollary 1 and 2.

A.1 Matrix concentration lemmas

We begin with matrices satisfying sub-Gaussianity of rows (Assumption 2) and then consider matrices whose rows are heavy-tailed (Assumption 1).

A.1.1 Supplemental lemmas for sub-Gaussian case

We define a sub-Gaussian random variable and a sub-Gaussian random vector373737Common examples of sub-Gaussian random variables are Gaussian, Bernoulli and all bounded random variables. below.

Definition 6**.**

A random variable $X$ is sub-Gaussian with parameter at most $K<\infty$ if for all $p\geq 1$, we have

[TABLE]

Further, a random vector $\mathbf{X}$ is sub-Gaussian with parameter at most $K$ if for every (fixed) vector $\mathbf{v}$, the random variable $\frac{\langle\mathbf{X},\,\mathbf{v}\rangle}{\|\mathbf{v}\|_{2}}$ is sub-Gaussian with parameter at most $K$.

We cite the following lemma for sub-Gaussian matrix concentration.

Lemma 3** (Theorem $5.39$,  [47]).**

Let $\mathbf{B}$ be a $n\times d$ matrix whose rows $\{\mathbf{b}_{i}\}_{i=1}^{n}$ are independent sub-Gaussian isotropic random vectors in $\mathbb{R}^{d}$. Then, for every $t>0$ we have the following inequality:

[TABLE]

with probability at least $1-2e^{-c_{K}t^{2}}$, where $C_{K},c_{K}$ depend only on the sub-Gaussian parameter $K$ of the columns.

To prove Corollary 1 for matrices $\mathbf{B}_{\mathsf{train}}$ satisfying the sub-Gaussian Assumption 2, we apply Lemma 3 for the matrix $\mathbf{B}:=\mathbf{B}_{\mathsf{train}}$ itself. We recall that we lower bounded the ideal test MSE as

[TABLE]

and then substituting Lemma 3 for the quantity $\sigma_{max}(\mathbf{B}_{\mathsf{train}})$ with $t:=\sqrt{n}$, together with the chi-squared tail bound on $\|\mathbf{W}_{\mathsf{train}}\|_{2}^{2}$ yields

[TABLE]

which is precisely Equation (4) with probability at least $(1-e^{-c_{K}n}-e^{-n\delta^{2}/8})$. This completes the proof of Corollary 1 for sub-Gaussian feature vectors. ∎

Before we move on to the heavy-tailed case, we remark on why we were not able to prove a version of Corollary 2 for the case of the feature vectors (rows of $\mathbf{B}_{\mathsf{train}}$) being sub-Gaussian. First, recall that to upper bound the ideal test MSE $\mathcal{E}_{\mathsf{test}}^{*}$, one needs to lower bound the minimum singular value of $\mathbf{B}_{\mathsf{train}}$ (or equivalently $\mathbf{B}_{\mathsf{train}}^{\top}$). Lemma 3 provides only a vacuous bound for the minimum singular value, as $n<d\leq C_{K}\sqrt{d}$ in general. Vershynin [47, Theorem $5.58$] does provide a concentration bound for tall matrices with independent, sub-Gaussian columns – which would be suitable to lower bound the minimum singular value of $\mathbf{B}_{\mathsf{train}}^{\top}$. However, this result uses an (unrealistically) restrictive condition of needing the columns to be exactly normalized as $\|\mathbf{b}_{i}\|_{2}=\sqrt{d}$ almost surely. Vershynin also provides an explicit example of a random tall matrix with sub-Gaussian columns that violates this condition, for which the minimum singular value does not concentrate. This suggests that sub-Gaussianity of high-dimensional feature vectors by itself may be too weak a condition to expect concentration on the minimum singular value, altogether a more delicate quantity.

However, we can prove Corollary 2 for the more special case where $\mathbf{B}_{\mathsf{train}}$ has independent, sub-Gaussian entries of unit variance. In other words, the random variable $B_{ij}$ is sub-Gaussian and unit variance, and the random variables $\{B_{ij}\}$ are independent. We cite the following lemma for concentration of the minimum singular value of such a random matrix $\mathbf{B}_{\mathsf{train}}$.

Lemma 4** (Theorem $5.38$, [47]).**

For $d\geq n$, let $\mathbf{B}$ be a $d\times n$ (or $n\times d$) random matrix whose entries are independent sub-Gaussian random variables with zero mean, unit variance, and sub-Gaussian parameter at most $K$. Then, for $\epsilon\geq 0$, we have

[TABLE]

where $C_{K}>0$ and $c_{K}\in(0,1)$ depend only on the sub-Gaussian parameter $K$.

We apply Lemma 4 substituting $\epsilon=\frac{1}{2C_{K}}$. Then, we get

[TABLE]

Noting as before that

[TABLE]

and substituting the upper chi-squared tail bound on the quantity $\|\mathbf{W}_{\mathsf{train}}\|_{2}^{2}$ as well as the lower tail bound on $\sigma_{min}(\mathbf{B}_{\mathsf{train}})$ above directly gives us the statement of Corollary 2 for the iid sub-Gaussian case.

A.1.2 Supplemental lemmas for heavy-tailed case

We cite the following lemma for heavy-tailed matrix concentration.

Lemma 5** (Theorem $5.41$,  [47]).**

Let $\mathbf{B}$ be a $n\times d$ matrix whose rows $\{\mathbf{b}_{i}\}_{i=1}^{n}$ are independent, isotropic random vectors in $\mathbb{R}^{d}$. Let $\|\mathbf{b}_{i}\|_{2}\leq\sqrt{n}$ almost surely for all $i\in[n]$. Then, for every $t>0$ we have the following inequality:

[TABLE]

with probability at least $1-2ne^{-ct^{2}}$ for some positive constant $c>0$.

We will use Lemma 5 to prove our lower bound for whitened feature matrices $\mathbf{B}_{\mathsf{train}}$ satisfying Assumption 1. Observe that Lemma 5 contains the deviation parameter $t$ as a multiplicative constant as opposed to an additive one, which is typical of concentration results for heavy-tailed distributions. Accordingly, the eventual scaling of the lower bound on the ideal test MSE, as well the extent of high probability in the result, will be accordingly weaker.

Substituting $t:=\sqrt{\frac{2}{c}(\ln n)}$ yields the upper tail inequality

[TABLE]

with probability at least $(1-\frac{1}{n})$. In a similar argument to the sub-Gaussian case, we can then lower bound the ideal test MSE as

[TABLE]

for constant $C>0$, completing the proof of Corollary 1 for heavy-tailed random feature vectors. ∎

As with the sub-Gaussian case, controlling the minimum singular value is much more delicate and in general it is not well-controlled. Vershynin [47, Theorem $5.65$] also has a result for tall matrices with independent columns more generally – but this result imposes even more conditions on top of the normalization requirement $\|\mathbf{b}_{i}\|_{2}=\sqrt{d}$: the columns need to be sufficiently incoherent in a pairwise sense.

A.2 Proof of Corollary 4

Recall from Section 5.1 that the basis pursuit program (BP) is given by

[TABLE]

Consider the following linear program (LP):

[TABLE]

Note that the concatenated vector $\begin{bmatrix}\mathbf{u}^{\top}&\mathbf{v}^{\top}\end{bmatrix}\in\mathbb{R}^{2d}$. Let $\widehat{\bm{\alpha}}$ refer to the optimal solution of BP in Equation (32) and let $\begin{bmatrix}\widehat{\mathbf{u}}^{\top}&\widehat{\mathbf{v}}^{\top}\end{bmatrix}$ be the optimal solution of the LP in Equation (33). A classical result is that the BP program in Equation (32) and the LP in Equation (33) are equivalent – not only in the sense that the maximal objective value is the same for both, but also that the optimal solution(s) are equivalent up to a linear transformation. We state the equivalence as a lemma and provide the entire proof argument for completeness.

Lemma 6**.**

The BP program (32) and the LP program (33) are equivalent in the following senses:

The maximal objectives are equal, i.e. $p^{*}_{BP}=p^{*}_{LP}$. 2. 2.

Given an optimal solution $(\widehat{\mathbf{u}},\widehat{\mathbf{v}})$ for the LP, we can obtain an optimal solution for the BP program through the linear transformation $\widehat{\bm{\alpha}}=\widehat{\mathbf{u}}-\widehat{\mathbf{v}}$. 3. 3.

Given an optimal solution $\widehat{\bm{\alpha}}$ for the the BP, we can obtain an optimal solution for the LP through the transformation

[TABLE]

Proof.

Denote the objective for the BP as $f(\bm{\alpha})=\|\bm{\alpha}\|_{1}$ and the objective for the LP as $g(\mathbf{u},\mathbf{v})=\sum_{i=1}^{d}\mathbf{u}_{i}+\sum_{i=1}^{d}\mathbf{v}_{i}$.

The crucial step is to show that a necessary condition for a solution $\begin{bmatrix}\mathbf{u}^{\top}&\mathbf{v}^{\top}\end{bmatrix}$ to be optimal is that $\mathsf{supp}(\mathbf{u})\cap\mathsf{supp}(\mathbf{v})=\emptyset$, i.e. that the supports of the vectors $\mathbf{u}$ and $\mathbf{v}$ are disjoint. This is stated formally in the lemma below.

Lemma 7**.**

The supports of any optimal solution $\begin{bmatrix}\widehat{\mathbf{u}}^{\top}&\widehat{\mathbf{v}}^{\top}\end{bmatrix}$ are disjoint, i.e $\mathsf{supp}(\widehat{\mathbf{u}})\cap\mathsf{supp}(\widehat{\mathbf{v}})=\emptyset$.

Taking Lemma 7 to be true for the moment, let us look at what it implies. For any solution where $\mathbf{u}$ and $\mathbf{v}$ are disjoint in their support, we can construct vector $\bm{\alpha}:=\mathbf{u}-\mathbf{v}$, and we note that

[TABLE]

In the other direction, for any feasible solution $\bm{\alpha}$, we can uniquely construct vectors $(\mathbf{u},\mathbf{v})$ as

[TABLE]

and we will again have

[TABLE]

Thus, we have established a one-one equivalence between all feasible solutions for the BP program and all feasible solutions for the LP that satisfy the disjoint-support condition. By Lemma 7, we know that all optimal solutions to the LP need to satisfy this condition. Thus, we have proved equivalence between the BP program (32) and the LP (33).

It only remains to prove Lemma 7, which we do below.

Proof.

It suffices to show that any solution not satisfying the disjointness property is strictly suboptimal. Suppose, for a candidate solution $\begin{bmatrix}\widehat{\mathbf{u}}^{\top}&\widehat{\mathbf{v}}^{\top}\end{bmatrix}$, there exists index $j\in[d]$ such that $\widehat{u}_{j}>0\text{ and }\widehat{v}_{j}>0$. Let $\widehat{c}_{j}=\widehat{u}_{j}-\widehat{v}_{j}$.

If $\widehat{c}_{j}>0$, the alternate solution $\begin{bmatrix}\mathbf{\widetilde{u}}^{\top}&\mathbf{\widetilde{v}}^{\top}\end{bmatrix}$ with $\widetilde{u}_{i}=\widehat{u}_{i},\widetilde{v}_{i}=\widehat{v}_{i}$ for $i\neq j$ and $\widetilde{u}_{j}=\widehat{c}_{j},\widetilde{v}_{j}=0$ will be feasible and have lower objective values. To see this, observe that

[TABLE]

where the last step follows since $\widehat{u}_{j},\widehat{v}_{j}>0$. On the other hand, if $\widehat{c}<0$, we can construct the alternative solution the alternate solution $\begin{bmatrix}\mathbf{\widetilde{u}}^{\top}&\mathbf{\widetilde{v}}^{\top}\end{bmatrix}$ with $\widetilde{u}_{i}=\widehat{u}_{i},\widetilde{v}_{i}=\widehat{v}_{i}$ for $i\neq j$ and $\widetilde{u}_{j}=0,\widetilde{v}_{j}=-\widehat{c}_{j}$. Clearly, this solution is feasible and by a similar argument also has lower objective value. This completes the proof. ∎

Since we have proved Lemma 7, we have completed our proof of equivalence. ∎

As a corollary of the above equivalence, it is clear that the support of any optimal solution for the BP program, $\widehat{\bm{\alpha}}$ is the same as the support of its equivalent solution $\begin{bmatrix}\widehat{\mathbf{u}}^{\top}&\widehat{\mathbf{v}}^{\top}\end{bmatrix}$ that is optimal for the LP. That is, we have

[TABLE]

As a consequence of the BP-LP equivalence, we can show through basic LP theory that there exists an optimal solution for the BP program whose support size is at most $n$. We state and prove this lemma below.

Lemma 8**.**

*There exists an optimal solution for the BP program (32) for which the support size is at most $n$.

Proof.

First we show that there exists an optimal solution $\begin{bmatrix}\widehat{\mathbf{u}}^{\top}&\widehat{\mathbf{v}}^{\top}\end{bmatrix}$ to the equivalent LP (33) with support size at most $n$. A basic feasible solution (BFS) for the LP is of the form $\begin{bmatrix}\mathbf{u}^{\top}&\mathbf{v}^{\top}\end{bmatrix}$ with support size at most $n$. It is well-known that for any LP, basic feasible solutions are extreme points of the feasible set [61]. Now, from Bauer’s maximum principle, we know that the minimum of a linear objective over a convex compact set is attained some extreme point of the feasible set, i.e. a BFS. Thus, there exists an extreme point, i.e. BFS, $\begin{bmatrix}\widehat{\mathbf{u}}^{\top}&\widehat{\mathbf{v}}^{\top}\end{bmatrix}$ that attains the minimum objective value for the LP and is thus an optimal solution. By the definition of BFS, this solution has support size at most $n$. Moreover, by the second statement of Lemma 6, $\widehat{\bm{\alpha}}=\widehat{\mathbf{u}}-\widehat{\mathbf{v}}$ is an optimal solution for the BP. From Equation (34) we have

[TABLE]

This completes the proof. ∎

Clearly, there exists an optimal solution of support size at most $n$. Generically, this optimal solution will be unique, and its support size will be exactly $n$. Showing these two properties is sufficient to show that BP is $1$-parsimonious (as defined in Definition 4). We conclude this section by proving these properties one-by-one for an appropriately non-degenerate matrix $\mathbf{A}_{\mathsf{train}}$, and Gaussian noise $\mathbf{W}_{\mathsf{train}}\sim\mathcal{N}(\mathbf{0},\sigma^{2}\mathbf{I}_{n})$. The sense in which we define non-degeneracy of covariate matrix $\mathbf{A}_{\mathsf{train}}$ is in the following three assumptions, which are extremely weak and will in general hold for any random ensemble with probability $1$. For any subset $S\subset[d]$, we denote the submatrix of $\mathbf{A}_{\mathsf{train}}$ corresponding to columns in subset $S$ by $\mathbf{A}_{\mathsf{train}}(S)$.

Assumption 4**.**

For every subset $S$ of size $n$, the matrix $\mathbf{A}_{\mathsf{train}}(S)\neq\mathbf{O}$.

Assumption 5**.**

For every subset $S$ of size $n$, the matrix $\mathbf{A}_{\mathsf{train}}(S)$ is invertible.

Assumption 6**.**

For any two distinct subsets $S_{1},S_{2}$ of size $n$, the matrix $(\mathbf{A}_{\mathsf{train}}(S_{1}))^{-1}\mathbf{A}_{\mathsf{train}}(S_{2})\notin\mathcal{P}^{n}$, where $\mathcal{P}^{n}$refers to the set of generalized permutation matrices of size $n\times n$.

Assumption 6 will be used to show uniqueness of the optimal solution to BP, which we state below as a lemma.

Lemma 9**.**

*The solution $\widehat{\bm{\alpha}}$ is unique with probability 1 over Gaussian noise $\mathbf{W}_{\mathsf{train}}\sim\mathcal{N}(0,\sigma^{2}\mathbf{I}_{n})$ when $\mathbf{A}_{\mathsf{train}}$ satisfies Assumptions 4, 5 and 6. *

Proof.

To show uniqueness of the optimal solution to the BP program (32), it suffices to show that there cannot exist more than one optimal BFS for the equivalent LP (33). We will show that the probability that there exists two basic feasible solutions for the LP that attain the same objective value is exactly zero. Let $\begin{bmatrix}\mathbf{u}_{1}^{\top}&\mathbf{v}_{1}^{\top}\end{bmatrix}$ and $\begin{bmatrix}\mathbf{u}_{2}^{\top}&\mathbf{v}_{2}^{\top}\end{bmatrix}$ be two basic feasible solutions to the LP with supports restricted to $S_{1}$ and $S_{2}$ where $S_{1},S_{2}\subset[d]$ and $S_{1}\neq S_{2}$. Let $\bm{\alpha}_{1}=\mathbf{u}_{1}-\mathbf{v}_{1}$ and $\bm{\alpha}_{2}=\mathbf{u}_{2}-\mathbf{v}_{2}$ denote their equivalent solutions for the BP program. Since $|S_{1}|=|S_{2}|=n$, from Assumption 5 we can write $\bm{\alpha}_{1}=(\mathbf{A}_{\mathsf{train}}(S_{1}))^{-1}\mathbf{W}_{\mathsf{train}},\bm{\alpha}_{2}=(\mathbf{A}_{\mathsf{train}}(S_{2}))^{-1}\mathbf{W}_{\mathsf{train}}.$ Denoting $\mathbf{V}:=(\mathbf{A}_{\mathsf{train}}(S_{2}))^{-1}\mathbf{W}_{\mathsf{train}}$, we have $\mathbf{V}\sim\mathcal{N}(0,(\mathbf{A}_{\mathsf{train}}(S_{2}))^{-1}((\mathbf{A}_{\mathsf{train}}(S_{2}))^{-1})\top)$. Thus, we get

[TABLE]

From the well-known result on the isometry of the $\ell_{1}$ norm [62], the last equality holds iff either $\mathbf{V}=\mathbf{0}$ or if $(\mathbf{A}_{\mathsf{train}}(S_{1}))^{-1}\mathbf{A}_{\mathsf{train}}(S_{2})\in\mathcal{P}^{n}$. From our assumptions,$(\mathbf{A}_{\mathsf{train}}(S_{1}))^{-1}\mathbf{A}_{\mathsf{train}}(S_{2})\notin\mathcal{P}^{n}$ and the equality can hold only if $\mathbf{V}=\mathbf{0}$. However, by Assumption 4, the random variable $\mathbf{V}$ has a non-zero covariance matrix, and so $\Pr\left[\|\bm{\alpha}_{1}\|_{1}=\|\bm{\alpha}_{2}\|_{1}\right]=0$. Hence, we have $\Pr\left[\|\bm{\alpha}_{1}\|_{1}=\|\bm{\alpha}_{2}\|_{1}\right]=0$ and the probability that two BFS’s can both be optimal for a random realization of noise $\mathbf{W}_{\mathsf{train}}$ is equal to [math]. This completes the proof. ∎

Finally, we show that the unique solution is of support size exactly $n$, which establishes $1$-parsimony for BP and allows us to apply Theorem 2 for the $\ell_{1}$-minimizing interpolator.

Lemma 10**.**

Under the above assumptions, the support of the unique solution $\widehat{\bm{\alpha}}$ is exactly $n$ with probability $1$ over Gaussian noise $\mathbf{W}_{\mathsf{train}}\sim\mathcal{N}(0,\sigma^{2}\mathbf{I}_{n})$.

Proof.

It is sufficient to show that for every $k<n$, $\Pr\left[|\mathsf{supp}(\widehat{\bm{\alpha}})|=k\right]=0$. Consider any subset $S\subset[d]$, $|S|=k$. Denote the columns of $\mathbf{A}_{\mathsf{train}}(S)$ as $\mathbf{a}_{1},\mathbf{a}_{2},\dots,\mathbf{a}_{k}$. Since $k<n$, there exists a vector $\mathbf{u}$ of unit $\ell_{2}$ norm such that $\mathbf{u}^{\top}\mathbf{a}_{i}=0$ for all $i\in[k]$. Let $\mathcal{E}_{S}$ denote the event that $\mathbf{W}_{\mathsf{train}}$ can be expressed as a linear combination of columns of $\mathbf{A}_{\mathsf{train}}(S)$. If $\mathcal{E}_{S}$ occurs then $\mathbf{W}_{\mathsf{train}}$ must be orthogonal to $\mathbf{u}$ i.e. $\mathbf{W}_{\mathsf{train}}^{\top}\mathbf{u}=0$. Thus,

[TABLE]

where the last equality follows by noting that $\|\mathbf{u}\|_{2}=1$ and thus $\mathbf{W}_{\mathsf{train}}^{\top}\mathbf{u}\sim\mathcal{N}(0,\sigma^{2})$. This argument clearly holds for any subset $S\subset[d],|S|=k$, and for any $k<n$. Therefore, the union bound gives us

[TABLE]

and since we already know that $|\mathsf{supp}(\widehat{\bm{\alpha}})|\leq n$, we have $\Pr\left[|\mathsf{supp}(\widehat{\bm{\alpha}})|=n\right]=1$, thus completing this proof. ∎

A.3 Proof of Proposition 1

In this section, we provide the proof of Proposition 1; in particular, the statement of Equation (29). Recall that we have assumed $\mathbb{E}\left[\mathbf{a}(\mathbf{X})\mathbf{a}(\mathbf{X})^{\top}\right]=\bm{\Sigma}=\mathbf{I}_{d}$. Using the triangle inequality, we have

[TABLE]

From the definition of the estimator $\widehat{\bm{\alpha}}_{\mathsf{hybrid}}$ and a similar argument as in the proof of Theorem 1, we have $\widehat{\bm{\alpha}}_{\mathsf{hybrid}}-\widehat{\bm{\alpha}}_{1}=\widehat{\bm{\Delta}}=\mathbf{A}_{\mathsf{train}}^{\dagger}\mathbf{W}_{\mathsf{train}}^{\prime}$ and thus,

[TABLE]

Recalling that we denoted $\mathcal{E}_{\mathsf{est}}(\widehat{\bm{\alpha}}_{1}):=\|\widehat{\bm{\alpha}}_{1}-\bm{\alpha}^{*}\|_{2}^{2}$, we have

[TABLE]

Now, remember that

[TABLE]

and thus we have

[TABLE]

where we recall the definition of the prediction error $\mathcal{E}_{\mathsf{est}}(\widehat{\bm{\alpha}}_{1})$. Substituting this in the above expression for the tesr error completes the proof of Equation (29) and thus the proof of Proposition 1.

∎

Recall that $\widehat{\bm{\alpha}}_{1}$ is a suitable estimator that uses the sparsity level $k$ or the noise variance $\sigma^{2}$ to estimate $\bm{\alpha}^{*}$ directly. So $\widehat{\bm{\alpha}}_{1}$ can be provided by a suitable estimator used for sparse recovery (like Lasso or OMP, but also many other choices). The results were stated as corollaries in Section 5.2. We provide the proofs of those corollaries below.

A.3.1 Proof of Corollary 5 (order-optimality)

Substituting the conditions for order-optimality (Equations (30a) and 30b) into Equation (29), we get

[TABLE]

In the proof of Corollary 5, we proved that $\|\mathbf{W}_{\mathsf{train}}\|_{2}^{2}\geq n\sigma^{2}(1-\delta)$ and $\lambda_{min}(\mathbf{A}_{\mathsf{train}}\mathbf{A}_{\mathsf{train}}^{\top})\geq(\sqrt{d}-2\sqrt{n})^{2}$ with probability at least $(1-e^{-n\delta^{2}/8}-e^{-n/2})$. Substituting these inequalities above, we get

[TABLE]

where in the inequality we used $d\geq cn$ for some $c>4$, which implies that $\sqrt{d}\geq\sqrt{c}\sqrt{n}\implies\sqrt{d}-2\sqrt{n}\geq(\sqrt{c}-2)\sqrt{n}=c^{\prime}\sqrt{n}$ for some $c^{\prime}>0$. This matches the statement of Equation (31) and completes the proof.

A.3.2 Recovery guarantees for Lasso (Corollary 6) and square-root-Lasso (Corollary 7)

In this section, we provide corollaries for hybrid interpolators that use the Lagrangian Lasso and the square-root Lasso. We follow notation from Wainwright’s book on high-dimensional non-asymptotic statistics [49, Chapter $7$] and provide original citations where-ever applicable.

Definition 7** (Lagrangian Lasso in presence of noise).**

We define the Lagrangian lasso with regularization parameter $\lambda_{n}$ (which can in general depend on $\sigma^{2}$ and $n$) as the estimator that solves the following optimization problem:

[TABLE]

Definition 8** (Square-root-Lasso).**

We define the square-root-Lasso with regularization parameter $\gamma_{n}$ (which can in general depend on $n$) as the estimator that solves the following optimization problem:

[TABLE]

For recovery guarantees on both the Lagrangian Lasso and the square-root-Lasso, in addition to $k$-sparsity, we require the following restricted eigenvalue condition on the design matrix:

Definition 9**.**

The matrix $\mathbf{A}_{\mathsf{train}}$ satisfies the restricted eigenvalue condition over set $\mathsf{supp}(\bm{\alpha}^{*})$ with parameters $(\kappa,\beta)$ if

[TABLE]

where for any subset $S$, we define

[TABLE]

The following lemma by Raskutti, Wainwright and Yu [63] shows that the matrix $\mathbf{A}_{\mathsf{train}}$ with iid standard normal entries (thus satisfying Assumption 3 atisfies this property with high probability as restated in the following lemma383838In fact the original result in the paper shows this property for generally correlated Gaussian design..

Lemma 11** (Corollary 1, [63]).**

The iid Gaussian matrix $\mathbf{A}_{\mathsf{train}}$ satisfies the restricted eigenvalue condition $(\kappa,\beta=3)$ for any $\kappa>0$ with probability greater than or equal to $(1-2e^{-cn})$, provided that $n\geq C\frac{k\ln d}{\kappa^{2}}$.

Under this condition, Bickel, Ritov and Tsybakov [56] proved the following bounds on estimation error as well as prediction error of the Lagrangian Lasso.

Theorem 3** (Theorems $7.13$ and $7.20$, [49]).**

Let $\mathbf{A}_{\mathsf{train}}$ satisfy the restricted eigenvalue condition with parameters $(\kappa,\beta=3)$. Then, any solution of the Lagrangian Lasso with regularization parameter $\lambda_{n}\geq 2\|\frac{\mathbf{A}_{\mathsf{train}}^{\top}\mathbf{W}_{\mathsf{train}}}{n}\|_{\infty}$ satisfies

[TABLE]

As a corollary, for any $\delta>0$ and regularization parameter choice $\lambda_{n}=2C^{\prime}\sigma\left(\sqrt{\frac{2\ln d}{n}}+\delta\right)$ we have

[TABLE]

with probability greater than or equal to $(1-2e^{-n\delta/2})$.

Observe from Equation (36a) and (36b) that $\mathcal{E}_{\mathsf{est}}\sim\mathcal{E}_{\mathsf{pred}}\leq C^{\prime\prime}\frac{\sigma^{2}k\ln d}{n}$ with high probability. Thus, a similar argument as in the proof of Corollary 5 will also give the statement of Corollary 6 for the case of the Lasso. We omit the argument for brevity. ∎

Similarly, the following guarantee was proved for the square-root-Lasso by Belloni, Chernozhukov and Wang [59], again assuming the restricted eigenvalue condition.

Theorem 4** ( [59]).**

Let $\mathbf{A}_{\mathsf{train}}$ satisfy the restricted eigenvalue condition with parameters $(\kappa,\beta=3)$. Then, any solution of the square-root-Lasso with regularization parameter $\gamma_{n}\geq 2\frac{\|\mathbf{A}_{\mathsf{train}}^{\top}\mathbf{W}_{\mathsf{train}}\|_{\infty}}{\sqrt{n}\|\mathbf{W}_{\mathsf{train}}\|_{2}}$ satisfies

[TABLE]

for constants $C,C^{\prime}>0$. As a corollary, for any $\delta>0$ and regularization parameter choice $\gamma_{n}=2C^{\prime}\left(\sqrt{\frac{2\ln d}{n}}+\delta\right)$ we have

[TABLE]

with probability greater than or equal to $(1-2e^{-n\delta/2})$ for new constants $C^{\prime\prime},C^{\prime\prime\prime}>0$.

Observe that the above choice of regularization parameter $\gamma_{n}$ does not depend on the noise variance $\sigma$ – thus, the square root Lasso can be implemented successfully, giving the rates in Equations (38a) and 38b, without requiring this knowledge. As in the case of the Lagrangian Lasso, we have $\mathcal{E}_{\mathsf{est}}\sim\mathcal{E}_{\mathsf{pred}}\leq C^{\prime\prime}\frac{\sigma^{2}k\ln d}{n}$ with high probability. Thus, a similar argument as in the proof of Corollary 5 will also give the statement of Corollary 7 for the case of the square-root-Lasso. This completes the proof of Corollary 7. ∎

A.3.3 Recovery guarantee for OMP (Corollary 6)

In this section we provide a formal statement for the recovery guarantee for orthogonal matching pursuit (OMP) in Corollary 6. The first theoretical guarantees on OMP were proved for noiseless recovery, when OMP is run to completion in the absence of noise [64]. Subsequent work [65] improved the sampling requirement, also in the noiseless setting. Recovery guarantees in the presence of noise were then proved for an appropriate stopping condition [57]. We formalize the version of OMP that we use below.

Definition 10** (OMP for recovery in presence of Gaussian noise [57]).**

Denote that $j^{th}$ column of $\mathbf{A}_{\mathsf{train}}$ by $\mathbf{a}_{j}$. The OMP algorithm is defined iteratively according to the following steps:

(Iteration $t=1$.) Initialize residual $r_{0}=\mathbf{Y}_{\mathsf{train}}$ and initialize the set of selected variables $S=\phi$. 2. 2.

Select variable $s_{t}:={\arg\max}|\langle\mathbf{a}_{j},\,\mathbf{r}_{t-1}\rangle|$ and add it to $S$. 3. 3.

Denote $\mathbf{A}_{\mathsf{train}}(S)$ as the submatrix of $\mathbf{A}_{\mathsf{train}}$ whose columns are in $S$. Let $\mathbf{P}_{t}=\mathbf{A}_{\mathsf{train}}(S)\mathbf{A}_{\mathsf{train}}(S)^{\dagger}$ denote the projection of $\mathbf{Y}_{\mathsf{train}}$ onto the linear space spanned by the elements of $\mathbf{A}_{\mathsf{train}}(S)$. Update $\mathbf{r}_{t}=(\mathbf{I}_{d}-\mathbf{P}_{t})\mathbf{Y}_{\mathsf{train}}$. 4. 4.

We stop the algorithm under one of two conditions: a) If $\|\mathbf{A}_{\mathsf{train}}\mathbf{r}_{t}\|_{\infty}\leq\sigma\sqrt{2(1+\eta)\ln d}$ for algorithmic parameter $\eta>0$; b) if the number of steps is equal to $k_{0}$, where $k_{0}>k$ is guaranteed. Otherwise, set $t=t+1$ and go back to Step 2. We refer to the stopping conditions henceforth as Condition a) and Condition b) respectively. 5. 5.

(Once algorithm has terminated) return $\widehat{\bm{\alpha}}_{1,\mathsf{OMP}}(\mathbf{A}_{\mathsf{train}},\mathbf{Y}_{\mathsf{train}})=\widehat{\bm{\alpha}}_{\mathsf{OLS}}(\mathbf{A}_{\mathsf{train}}(S),\mathbf{Y}_{\mathsf{train}})$.

We also define pairwise incoherence

[TABLE]

In general for recovery guarantees, we would like this quantity to be small. We state the following well-known lemma for the entries of $\mathbf{A}_{\mathsf{train}}$ being iid $\mathcal{N}(0,1)$:

Lemma 12**.**

For iid Gaussian $n\times d$ matrix $\mathbf{A}_{\mathsf{train}}$, we have $\mu(\mathbf{A}_{\mathsf{train}})<\frac{1}{2k-1}$ with high probability as long as $n=\Omega(k^{2}\ln d)$.

We now restate the result that is most relevant for our purposes, and holds for any design matrix $\mathbf{A}_{\mathsf{train}}$.

Theorem 5** (Theorem $8$ from [57].).**

Suppose that $\mu<\frac{1}{2k-1}$ and all the non-zero coefficients $\alpha^{*}_{i},i\in\mathsf{supp}(\bm{\alpha}^{*})$ satisfy

[TABLE]

for some $\eta>0$. Then, the above version of OMP has the following guarantees:

Under stopping condition a), OMP will recover $S=\mathsf{supp}(\bm{\alpha}^{*})$ with probability greater than or equal to $(1-\frac{k}{d^{\eta}\sqrt{2\ln d}})$. This directly implies that

[TABLE]

with probability greater than or equal to $\left(1-\frac{k}{d^{\eta}\sqrt{2\ln d}}\right)$. 2. 2.

Under stopping condition b), OMP will recover $S\supset\mathsf{supp}(\bm{\alpha}^{*})$ with probability greater than or equal to $(1-\frac{k_{0}}{d^{\eta}\sqrt{2\ln d}})$. This directly implies that

[TABLE]

with probability greater than or equal to $\left(1-\frac{k_{0}}{d^{\eta}\sqrt{2\ln d}}\right)$.

We observe that the obtained guarantees for OMP in the noisy setting are support recovery guarantees, which are much stronger than directly bounding the $\ell_{2}$-norm of the error. While the expressions do not have an implicit dependence on the dimension $d$, it requires the absolute values of the coefficients $|\alpha^{*}_{i}|$ to be bounded above by $\Omega(\sigma\sqrt{\ln d})$. This can effectively be thought of as requiring the SNR to be high enough to ensure signal recovery - the more we overparameterized the number of features $d$, the more stringent this condition becomes. Also observe that the versions of OMP with stopping conditions a) and b) use different kinds of side information: Stopping condition a) uses knowledge of the noise variance, while stopping condition b) only uses a guaranteed upper bound on the true sparsity level $k$. Stopping condition a) gives the bound in Equation (40a) which is adaptive to unknown sparsity level. Stopping condition b) has a weaker kind of side information, but gives the bound in Equation (41a) which, while having a good convergence rate, does not adapt to the unknown sparsity level.

As before, a similar argument as in the proof of Corollary 5 will also give the statement of Corollary 6 for the case of OMP with both stopping conditions. We omit the argument for brevity.

Appendix B The interpolation threshold: Regularity vs randomness of training data

We observe that Corollary 2 is meaningful primarily for the heavily overparameterized regime where $d>>n$ (more formally, if we vary $d$ as a function of $n$, we have $\lim_{n\to\infty}\frac{n}{d(n)}=0$). It does not mathematically explain the magnitude of the interpolation peak at $d\sim n$ that we observe in Figure 1, which reflects the mean of the minimum possible test MSE that results purely from fitting noise at the interpolation threshold. It is well-known that the behavior of the minimum singular value of a random matrix with iid Gaussian entries is very different when the matrix is approximately square – a “hard-edge" phenomenon arises [66, 67] and the minimum singular value actually exhibits heavy-tailed behavior [68]. It is unclear whether the same phenomena hold for lifted features on data such as Fourier, Legendre or even Vandermode features. It thus makes sense to consider the statistics of the test MSE at the interpolation threshold more carefully.

Figure 9 shows the generalization performance of various interpolating solutions on approximately “whitened" Legendre and Fourier features on data. In particular, we consider two generative assumptions for the training data $\{X_{i}\}_{i=1}^{n}$:

Regularly spaced data: $X_{i}=-1+\frac{2(i-1)}{(n-1)}\text{ for all }i\in[n]$. 2. 2.

Randomly drawn data: $X_{i}\sim\text{Unif}[-1,1]$.

and the generative model for the data to be $Y_{i}=f^{*}(X_{i})+W_{i}$ where $f^{*}:[-1,1]\to\mathbb{R}$ is a polynomial of degree $2$. We make $50$ random draws of both the training data (in the case of randomly drawn training points) and noise, and plot both the mean and the median test MSE evaluated over $10^{4}$ test points, together with $85\%$ confidence intervals. The results depicted in Figure 9 have several implications. In particular, the distributional assumptions on the data make a significant difference: regularly spaced data approximately preserves the orthogonality of the design matrix $\mathbf{A}_{\mathsf{train}}\approx\mathbf{B}_{\mathsf{train}}$393939recall that the Legendre and Fourier features are orthogonal in function space, which corresponds to an integral over the interval $[-1,1]$; regular spacing of training points appropriately approximates this integral.. On the other hand, randomly drawn training points constitute a more complex random matrix $\mathbf{A}_{\mathsf{train}}$, for which the minimum singular value is much lower and moreover exhibits heavy-tailed behavior404040Intuitively, this is because of connections of the minimum singular value of the random matrix to the non-uniformity of random spacings [68]; one can prove that that the smallest spacing between randomly drawn points is on the order of $1/n^{2}$ (as opposed to $1/n$) with high probability. . This can be observed from the large size of the confidence intervals, as well as the much lower median test MSE (which is still higher than the test MSE incurred from regularly spaced training data). Randomly drawn training data is what we use in the conventional machine learning step, and this analysis tells us that we truly need to overparameterize to ensure harmless fitting of noise, unlike in logistic regression.

Appendix C Mathematical facts

In this section, we collect miscellaneous mathematical facts that were useful for proving some of our results.

Fact 1**.**

For positive numbers $(n,k)$, we have

[TABLE]

Proof.

Observe that

[TABLE]

Next, we have $e^{k}=\sum_{i=0}^{\infty}\frac{k^{i}}{i!}\geq\frac{k^{k}}{k!}$. Rearranging this gives us $k!\geq\frac{k^{k}}{e^{k}}$, and substituting it above gives us

[TABLE]

completing the proof of this fact. ∎

Appendix D Calculations for the regularly spaced Fourier case

Consider observation $\mathbf{y}_{\mathsf{train}}\in\mathbb{C}^{n}$ and $d$ Fourier features, $f_{k}(x_{\mathsf{train,j}})=e^{i2\pi kx_{\mathsf{train,j}}}$ for $k=\{0,1,\dots,d-1\}.$ Here $\mathbf{x}_{\mathsf{train}}$ contains regularly spaced samples with $x_{\mathsf{train,j}}=\frac{j}{n}$ for $j\in\{0,1,\dots,n-1\}.$ For ease of exposition we assume $d=(M+1)n$ for positive integer $M$. We wish to understand how the minimum weighted $\ell_{2}$ norm interpolating solution behaves in this setting. The first $n$ features, form an orthogonal basis for $\mathbb{C}^{n}$ and thus it suffices to understand the case where $\mathbf{y}_{\mathsf{train}}$ is each basis vector separately. Let

[TABLE]

for some $k^{*}\in\{0,1,\dots,n-1\}$.Without loss of generality we consider $k^{*}$ in the range $[0,n-1]$ since subsequent blocks of $n$ features will be aliases of these features.

Note that we can write,

[TABLE]

where $\boldsymbol{\alpha^{*}}=\mathbf{e_{k^{*}}}\in\mathbb{R}^{n}$ and $\mathbf{e}_{j}$ denotes the $j^{th}$ standard basis vector. For any solution $\boldsymbol{\alpha}\in\mathbb{C}^{n}$, the interpolating constraint is,

[TABLE]

If we scale feature $f_{j}(\mathbf{x}_{\mathsf{train}})$ by real weight $w_{j}$, then the interpolating constraint becomes,

[TABLE]

with $\beta_{j}=\frac{\alpha_{j}}{w_{j}}$.

We are interested in the minimum weighted $\ell_{2}$ norm solution subject to the interpolating constraint given by,

[TABLE]

where $\boldsymbol{\Gamma}=\text{diag}(w_{0},w_{1},\dots,w_{d-1})$. Note that this is equivalent to the minimum $\ell_{2}$-minimizing coefficients corresponding to the weighted featurues, as defined in [8]:

[TABLE]

We will solve the problem in (44) next. First, we list some properties of the regularly spaced Fourier features. Denote by $S(k^{*})$, the set of indices corresponding to features that are exact aliases of $f_{k^{*}}(\mathbf{x}_{\mathsf{train}})$. Then,

[TABLE]

Note that $M=|S(k^{*})|=\frac{d}{n}-1$. We have,

[TABLE]

Using (42), (46) and (47) we can rewrite the optimization problem in (44) as,

[TABLE]

Clearly to minimize the objective we must have $\hat{{\beta}}_{j}=0$ for $j\notin\{\{k^{*}\}\cup S(k^{*})\}$ and thus it suffices to consider the problem restricted to the indices in $\{\{k^{*}\}\cup S(k^{*})\}$. By mapping these indices to the set $\{0,1,\dots,M\}$ and denoting the weight vector restricted to this set as $\mathbf{\tilde{w}}$ we write an equivalent optimization problem,

[TABLE]

Note that $\tilde{w}_{j}=w_{k^{*}+jn}$ for $j=0,1,\dots,M$.

We find an optimal solution to this problem by using the Cauchy Schwarz inequality which states,

[TABLE]

where equality occurs if and only if $\hat{\boldsymbol{\xi}}=c\mathbf{\tilde{w}}$ for some $c\in\mathbb{C}$. Using the fact that $w_{j}\in\mathbb{R}$ and solving for $c$ using the interpolating constraint we get,

[TABLE]

Mapping this back to original indices we get,

[TABLE]

and

[TABLE]

where,

[TABLE]

Next we consider the effect on a test point $\mathbf{X}\in\mathbb{R}$ with i.i.d. entries $\mathbf{X}\sim U[0,1]$, when the ground truth observation is $\mathbf{Y}=f_{k^{*}}(\mathbf{X})$. On this point we predict,

[TABLE]

We want to understand how different $\mathbf{\hat{Y}}$ is from $\mathbf{Y}$. Using (48) we have,

[TABLE]

The prediction $\mathbf{\hat{Y}}$ consists of two components. The first component is the true signal attenuated by a factor $\hat{{\alpha}}_{k^{*}}$ due to the effect of signal bleed. The signal bleeds to features that are orthogonal to the true signal and this leads to the second component, a contamination term that we denote by,

[TABLE]

Let $\mathsf{SU}(k^{*})$ denote the fraction of the true coefficient that survives post signal bleed. Then,

[TABLE]

Let $C(k^{*})$ denote the standard deviation of the contamination given by,

[TABLE]

Using the property of Fourier features when $\mathbf{X}$ is spaced uniformly in $[0,1]$ namely,

[TABLE]

to get,

[TABLE]

Using this we have,

[TABLE]

Next we consider examples of weighting schemes for a given $n,d$ pair with large enough $\frac{d}{n}$ when the true signal is at $k^{*}$. The set of indices containing aliases of the true signal is denoted as $S(k^{*})$ as in (45).

Example 7**.**

Uniform weights, $w_{j}=1$.

[TABLE]

Example 8**.**

Spiked weights on low frequency features: This selects a fraction of energy to put on the favored set of $s<n$ features. Namely. for $\gamma\in[0,1]$ and $s<n$.

[TABLE]

For $0\leq k^{*}<s$,

[TABLE]

For $s\leq k^{*}<n$,

[TABLE]