Harmless interpolation of noisy data in regression

The ideal interpolator is the whitened minimum-norm fit of the noise (Theorem 1). With $\Delta=\alpha-\alpha^*$, the test error is $\|\Sigma^{1/2}\Delta\|^2$ and interpolation reads $A\Delta=W$. In whitened coordinates $u=\Sigma^{1/2}\Delta$ this is the minimum-norm point of $\{u:Bu=W\}$: $$u^\star=B^\top(BB^\top)^{-1}W,\qquad \|u^\star\|^2=W^\top(BB^\top)^{-1}W.$$ Nothing random is used. Each extra feature adds a PSD term to $BB^\top$ and lowers this minimum; with $d\gg n$, $BB^\top\approx dI_n$ and the noise energy $n\sigma^2$ is divided among $d$ weak directions.

The price of interpolation is bounded by the extreme singular values of $B$, and the lower bound on the smallest one is the hard part (Corollaries 1 and 2). Both corollaries come from $\|W\|^2/\sigma_{\max}(B)^2\le W^\top(BB^\top)^{-1}W\le\|W\|^2/\sigma_{\min}(B)^2$. Bounding $\sigma_{\max}$ above is easy (Davidson-Szarek, Vershynin, and matrix Bernstein for bounded rows, which introduces the $\sqrt{\ln d}$). Bounding $\sigma_{\min}$ below is delicate: the Gaussian tail gives $\sqrt d-2\sqrt n$ (hence $d>4n$) and Rudelson-Vershynin covers independent entries, but sub-Gaussian rows alone do not force $\sigma_{\min}$ to concentrate, and for lifted features on random points the paper's experiments find it heavy-tailed near $d\approx n$.

On a regular grid, aliases absorb noise at exactly the ideal rate, and a true frequency survives in proportion to its share of the cohort's weight. Since $e^{2\pi i(k+\ell n)j/n}=e^{2\pi ikj/n}$, the $d$ features fall into $n$ cohorts of identical training columns, distinct cohorts are orthogonal on the sample, and all features are orthonormal at test time. Expand the noise as $W=\sum_{k<n}c_kf_k(x_{\rm train})$, $\mathbb{E}|c_k|^2=\sigma^2/n$. Minimum norm splits $c_k$ evenly over its $d/n$ cohort members (equality in Cauchy-Schwarz), so the total is $\|W\|^2/d\approx\sigma^2n/d$. Here $BB^*=dI_n$, so on pure noise the minimum-$\ell_2$ interpolator is the ideal one. With weights, the fit within a cohort is a matched filter: in rescaled coefficients $\xi_j=\alpha_j/w_j$ the problem is $\min\|\xi\|^2$ subject to $\langle\xi,w\rangle=1$, so $\xi=w/\|w\|^2$, $\hat\alpha_j=w_j^2/V$, and $\mathrm{SU}$ as above. The signal survives only if its weight dominates the total alias weight, which the weighting fixes before any data are seen; with uniform weights $\mathrm{SU}=n/d\to0$, which is why in the isotropic Gaussian case the minimum-$\ell_2$ interpolator tends to the null risk.

The energy bled into the aliases is the contamination. Exact interpolation puts the energy taken from the true feature on its aliases, which are orthogonal at test time and so add error there: orthonormality gives $\mathbb{E}|\hat Y-f_{k^*}(X)|^2=(1-\mathrm{SU})^2+C^2$ with $C^2=\sum_{S}\hat\alpha_j^2$. Contamination is a separate condition from bleed: uniform weights bleed completely yet have $C\approx\sqrt{n/d}\to0$, while a few heavy aliases keep $\mathrm{SU}$ away from 1 and $C$ away from 0.

Summing over cohorts gives the consistency conditions. Cohorts are orthogonal at training and test time and the fit is linear in $Y$, so for $\alpha^*$ supported on $\{0,\dots,n-1\}$, in expectation over the noise, $$\mathbb{E}\,\mathcal{E}_{\rm test}=\sum_{k^*<n}|\alpha^*_{k^*}|^2\big[(1-\mathrm{SU})^2+C^2\big]+\frac{\sigma^2}{n}\sum_{k^*<n}\big[\mathrm{SU}^2+C^2\big].$$ This assembly is mine; the paper states the conclusions in words. Consistency needs $\mathrm{SU}\to1$ and $C\to0$ where the signal lives and $\mathrm{SU}^2+C^2\to0$ on all but $o(n)$ cohorts, so the surviving frequencies must grow sub-linearly in $n$ and the rest must bleed out.

Spiked weights meet both conditions with a few strong features and many weak ones. For $k^*<s$ every alias is unfavored, and with $M\approx d/n$, $$\mathrm{SU}\approx\frac{1}{1+\frac sn\big(\frac1\gamma-1\big)},\qquad C\approx\frac{s/\sqrt{nd}}{\frac{\gamma}{1-\gamma}+\frac sn},$$ while for $k^*\ge s$ the cohort is uniform and $\mathrm{SU}^2+C^2=n/d$. The noise term is then about $\sigma^2(s/n+n/d)$ and the signal term, for a signal in the favored set, $((s/n)(1/\gamma-1))^2+O(s^2/(nd))$, so for fixed $\gamma\in(0,1)$ consistency holds iff $s/n\to0$ and $d/n\to\infty$ (my assembly; the paper gives only the per-cohort $\mathrm{SU}$ and $C$). The strong features carry the signal and also fit their own share of the noise; the weak ones absorb the rest.

Ridge is minimum-norm interpolation with $n$ extra features, so it reduces contamination but increases bleed. Eliminating the residual shows that ridge with penalty $\lambda^2\alpha^\top\Gamma\alpha$ is the minimum-norm interpolator on the augmented design $[A\Gamma^{-1/2},\ \lambda I_n]$. The ridge columns are indicators of the training points and vanish at test time, so energy bled into them causes no contamination; but they are aliases that need smaller coefficients as $\lambda$ grows, so they increase bleed. So optimally tuned ridge beats the minimum-$\ell_2$ interpolator (Hastie et al.) by reducing contamination, and only a reweighting $\Gamma$ aligned with the signal reduces bleed.

A parsimonious interpolator fits the noise with $n$ columns, and no $n$ columns have a large top eigenvalue (Theorem 2). Take $\alpha^*=0$, so $\mathcal{E}_{\rm test}=\|\hat\alpha\|^2\ge\|\hat\alpha_S\|^2$ on the top-$n$ support $S$. Parsimony gives $\|A_S\hat\alpha_S\|^2\ge\beta\|W\|^2$, and $\|A_S\hat\alpha_S\|^2\le\lambda_{\max}(A_S^\top A_S)\|\hat\alpha_S\|^2$, so $$\|\hat\alpha_S\|^2\ \ge\ \frac{\beta\|W\|^2}{\max_{|S'|=n}\lambda_{\max}(A_{S'}^\top A_{S'})}.$$ The maximum over all $n$-subsets is needed because $S$ depends on $W$. A Gaussian tail per block and a union bound over $\binom dn\le(ed/n)^n$ subsets give $\max\lambda_{\max}\approx4n\ln(d/n)$. In Step 1 all $d$ columns carry the noise and $\lambda\approx d$; here only $n$ do, and the logarithm from the union bound accounts for the whole $1/\ln(d/n)$ rate.

OMP to completion and basis pursuit are hard-sparse, so $\beta=1$. Both return interpolators on at most $n$ columns: OMP reaches zero residual after $n$ steps, and basis pursuit as an LP in $\alpha=u-v$ has a basic optimum with at most $n$ nonzeros, one per equality constraint. Uniqueness (or a solver that returns a vertex) makes the output this sparse optimum, and it holds a.s. when $A$ has a density, since a tie between optimal vertices forces an algebraic identity in $A$; the paper's extra claim of support exactly $n$ is not needed.

The hybrid's error splits orthogonally into missed signal and refitted noise (Proposition 1). The second stage is Step 1 applied to the residual. With $P$ the projection onto the row space of $A$, $$\hat\alpha_{\rm hyb}-\alpha^*=(I-P)(\hat\alpha_1-\alpha^*)+A^\top(AA^\top)^{-1}W,$$ and the two terms are orthogonal, so $\mathcal{E}_{\rm test}(\hat\alpha_{\rm hyb})\le\mathcal{E}_{\rm est}+\|W\|^2/\lambda_{\min}(AA^\top)$. The refit restores the row-space part of missed signal, so the signal term is at most $\mathcal{E}_{\rm est}$.

Plug in first-stage rates and match a lower bound (Corollaries 5 to 7). The hybrid's error is at most the first-stage rate plus the price of interpolation: SLOPE supplies $k\ln(d/k)/n$, the Lasso (under a restricted eigenvalue condition) and the $\sigma$-free square-root Lasso $k\ln d/n$, and OMP with a noise-level stopping rule $k\sigma^2/n$ under incoherence and a minimum signal size. Every interpolator has error at least the bound of Corollary 1, and every estimator at least the minimax rate on some $k$-sparse instance, and the max of the two is at least half their sum. As $d$ grows the noise term falls like $n/d$ while the signal term grows like $\ln d$, so the best $d$ is large but finite.