Harmless interpolation of noisy data in regression

Introduction

This wisdom has been challenged by the recent advent of deeper and deeper neural networks. In particular, a thought-provoking paper 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?

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 as the number of features goes to infinity for feature families satisfying mild conditions.

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).

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.

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 and ill-conditioning in least-squares regression . 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 . 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 , 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 reasonsThis 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 , is challenging.

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.

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 ) 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.

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 methodsfor which, interestingly, the induced inductive bias is unknown. in optimization need not always improve generalization for specially constructed examples in overparameterized linear regression . On the other hand, on a different set of examples , adaptive methods converge to better-generalizing solutions than SGD. Evidence suggests that norm-based complexity measures predict generalizing ability , and for neural networks, such complexity measures have been developed that do not depend on the width, but can depend on the depth . 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 and is recently given as a justification for using overparameterization in neural networks . However, margin does not always imply generalization . An alternative, intriguing explanation 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.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 empirically recovered several of the overparameterization phenomena of deep learning in kernel classifiers that interpolateIn 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 improvementThis was also observed by 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 , but most relevant to the setting of regression is recent analysis of the test MSE of the minimum-RKHS-norm kernel interpolator . 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 . 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.3 High-dimensional linear regression

1.4 Concurrent work in high-dimensional linear regression

Problem Setting

We define shorthand notation for the training data: let

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:

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

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.

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.

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

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.

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

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.

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

$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.

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

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

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

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

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

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.

[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.

[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.

[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.

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

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.

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$.

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 generalityNote 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$).

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.

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

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

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

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$.

$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 as $n\to\infty$ at a polynomial rate.

$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 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 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 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

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:

Observe that Equation (1) can be rewritten as

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}$:

(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:

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 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 . The original argument is contained in classical work by Davidson and Szarek .

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

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

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

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

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

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

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 variablesOne could have also used the Hanson-Wright inequality 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

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. ∎

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}^{*})=$. We take $n=5000$ measurements, each of which is corrupted by Gaussian noise of variance $\sigma^{2}=0.01$.

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:

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.

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.

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

The estimator has to interpolate this data with some linear combination of Fourier featuresWhy 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 aliasesThese 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$.

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 essentially 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}$.

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 priorThis 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 consideredIn 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.

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:

The output energy (signal as well as noise) bleeds away from the true signal component corresponding to Fourier feature – 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

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

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.

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 . 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.

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-frequencyThis 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.”

On one hand, if $\omega(n)$ of the $\mathsf{SU}$s stay boundedly above , 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.

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

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.

Interpolation in the noisy sparse linear model

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

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.

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

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, 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 as $d\to\infty$. This regime is plotted in Figure 8(a) for $n=1000$, and 8(b) for $n=50$.

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

$d=e^{n}$ and $n$ growing. As before, we want the test MSE to to decay to 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.

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

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

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).

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

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}}$.

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

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 as $n\to\infty$.

$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 as $n\to\infty$, but at an extremely slow logarithmic rate.

$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.

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

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.

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

We close this section with the proof of Theorem 2.

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

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

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.

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

Assumption 3 (Gaussian features), we have

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

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$,

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

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.

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

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

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

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. ∎

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 instancesStatisticans 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.

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.

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 expositionThis 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. .

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}$.

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.

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

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

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.

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

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

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

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 . 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.

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 ruleDetails 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

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

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.

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

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 discoveredThis 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 ., 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.

Conclusions for high-dimensional generative model

Our results, together with other work in the last year , 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 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.

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.

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 , or by first estimating the SNR .

Finally, as also pointed out by , 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 . 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.

References

Appendix A Supplemental proofs

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

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

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

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

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.

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

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

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}}$.

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

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

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.

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

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

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

Consider the following linear program (LP):

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}$.

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}}$.

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

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.

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

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

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.

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

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

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.

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

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 . 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

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)$.

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

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

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.

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.

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

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})$.

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

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. ∎