On Optimal Interpolation in Linear Regression
Year: 2021

Abstract:
Understanding when and why interpolating methods generalize well has recently been a topic of interest in statistical learning theory. However, systematically connecting interpolating methods to achievable notions of optimality has only received partial attention. In this paper, we investigate the question of what is the optimal way to interpolate in linear regression using functions that are linear in the response variable (as the case for the Bayes optimal estimator in ridge regression) and depend on the data, the population covariance of the data, the signal-to-noise ratio and the covariance of the prior for the signal, but do not depend on the value of the signal itself nor the noise vector in the training data. We provide a closedform expression for the interpolator that achieves this notion of optimality and show that it can be derived as the limit of preconditioned gradient descent with a specific initialization. We identify a regime where the minimum-norm interpolator provably generalizes arbitrarily worse than the optimal response-linear achievable interpolator that we introduce, and validate with numerical experiments that the notion of optimality we consider can be achieved by interpolating methods that only use the training data as input in the case of an isotropic prior. Finally, we extend the notion of optimal response-linear interpolation to random features regression under a linear data-generating model that has been previously studied in the literature.

heading: 1 Introduction
text: Establishing mathematical understanding for the good generalization properties of interpolating methods, i.e. methods that fit the training data perfectly, has attracted significant interest in recent years. Motivated by the quest to explain the generalization performance of neural networks which have zero training error, for example even on randomly corrupted data (Zhang et al., 2017), this area of research has established results for a variety of models. For instance, in kernel regression, Liang and Rakhlin (2020) provide a data-dependent upper bound on the generalization performance of the minimum-norm interpolator. By analyzing the upper bound, they show that small generalization error of the minimum-norm interpolator occurs in a regime with favourable curvature of the kernel, particular decay of the eigenvalues of the kernel and data population covariance matrices and, importantly, in an overparametrized setting. In random features regression, Mei and Montanari (2019) showed that for large signal-to-noise ratio and in the limit of large overparametrization, the optimal regularization is zero, i.e. the optimal ridge regressor is an interpolator. Liang and Sur (2020) characterized the precise high-dimensional asymptotic generalization of interpolating minimum`1-norm classifiers and boosting algorithms which maximize the `1 margin. Bartlett et al. (2020) isolated a setting of benign overfitting in linear regression, dependent on notions of effective rank of the population covariance matrix, in which the minimum-norm interpolator has small generalization error. Similarly, this regime of benign overfitting occurs with large overparametrization.
Linear models, in particular, provide a fundamental playground to understand interpolators. On the one hand, in overparametrized regimes, interpolators in linear models are seen to reproduce
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
stylized phenomena observed in more general models. For example, the double descent phenomenon, which was first empirically observed in neural networks (Belkin et al., 2019), has also featured in linear regression (Hastie et al., 2019). On the other hand, neural networks are known to be well-approximated by linear models in some regimes. For example, with specific initialization and sufficient overparametrization, two-layer neural networks trained with gradient descent methods are well-approximated by a first-order Taylor expansion around their initialization (Chizat et al., 2019). This linear approximation can be split into a random features component and a neural-tangent component. The random features model, a two-layer neural network with randomly initialized first layer which is fixed during training, shares similar generalization behavior with the full neural network (Bartlett et al., 2021), and as such, the random features model provides a natural stepping stone towards tackling a theoretical understanding of neural networks.
A major focus of the interpolation literature has so far been to theoretically study if and when interpolating methods based on classical techniques such as ridge regression and gradient descent can have optimal or near-optimal generalization (Bartlett et al., 2021). However the question of understanding which interpolators are best, and designing data-dependent schemes to implement them, seems to have received only partial attention. Work investigating which interpolators are optimal in linear regression includes Muthukumar et al. (2019), where the authors constructed the best-possible interpolator, i.e. a theoretical device which uses knowledge of the true parameter and training noise vector to establish a fundamental limit on how well any interpolator in linear regression can generalize. When the whitened features are sub-Gaussian, this fundamental limit is lower bounded by a term proportional to n/d, up to an additive constant and with high probability, which is small only in the regime of large overparametrization. Here, n and d are the size and the dimension of the data. While this interpolator provides the best-possible generalization error, the interpolator is not implementable in general, as one would need access to the realization of the true data-generating parameter w? and the realization of the noise in the training data. Rangamani et al. (2020) studied generalization of interpolators in linear regression and showed that the minimum-norm interpolator minimizes an upper bound on the generalization error related to stability. In (Mourtada, 2020), it was shown that the minimum-norm interpolator is minimax optimal over any choice of the true parameter w? ∈ Rd, distributions of the noise with mean 0 and bounded variance, and for a fixed nondegenerate distribution of the features. Amari et al. (2021) computed the asymptotic risk of interpolating preconditioned gradient descent in linear regression and investigated the role of its implicit bias on generalization. In particular, they identified the preconditioning which leads to optimal asymptotic (as d/n → γ > 1 with n, d → ∞) bias and variance, separately, among interpolators of the form w = PXT (XPXT )−1y for some matrix P , where X ∈ Rn×d is the data matrix, y ∈ Rn is the response vector. They showed that, within this class of interpolators, using the inverse of the population covariance matrix of the data as preconditioning achieves optimal asymptotic variance. However, the interpolator with optimal risk is not given.
In this paper, we study the question of what is the optimal way to interpolate in overparametrized linear regression by procedures that do not use the realization of the true parameter generating the data, nor the realization of the training noise. The motivation for studying this question is twofold. First, in designing new ways to interpolate that are directly related to notions of optimality in linear models, we hope to provide a stepping stone to designing new ways to interpolate in more complex models, such as neural networks. Second, our results illustrate that there can be arbitrarily large differences in the generalization performance of interpolators, in particular considering the minimum-norm interpolator as a benchmark. This is a phenomenon that does not seem to have received close attention in the literature and may spark new interest in designing interpolators connected to optimality.
We consider the family of interpolators that can be achieved as an arbitrary function f of the data, population covariance, signal-to-noise ratio and the prior covariance such that f is linear in the response variable y (as the case for the Bayes optimal estimator in ridge regression). We call such interpolators response-linear achievable (see Definition 3). We also introduce a natural notion of optimality that assumes that the realization of true data-generating parameter and the realization of the noise in the training data are unknown. Within this class of interpolators and under this notion of optimality, we theoretically compute the optimal interpolator and show that this interpolator is achieved as the implicit bias of a preconditioned gradient descent with proper initialization. We refer to this interpolator as the optimal response-linear achievable interpolator.
Could it be that the commonly used minimum-norm interpolator is good enough so that the benefit of finding a better interpolator is negligible? We illustrate that the answer to this question is, in general,
no. In particular, we construct an example in linear regression where the minimum-norm interpolator has arbitrarily worse generalization than the optimal response-linear achievable interpolator. Here, the variance (hence also generalization error) of the minimum-norm interpolator diverges to infinity as a function of the eigenvalues of the population covariance matrix, while the generalization error of the optimal response-linear achievable interpolator stays bounded, close to the optimal interpolator, i.e. the theoretical device of Muthukumar et al. (2019) which uses the value of the signal and noise.
The optimal response-linear achievable interpolator uses knowledge of the population covariance matrix of the data (similarly as in Amari et al. (2021)), the signal-to-noise ratio, and the covariance of the true parameter (on which we place a prior distribution). Is it the case that the better performance of our interpolator is simply a consequence of this population knowledge? We provide numerical evidence that shows that the answer to this question is, in general, no. In particular, we construct an algorithm to approximate the optimal response-linear achievable interpolator which does not require any prior knowledge of the population covariance or the signal-to-noise ratio and uses only the training data X and y, and we empirically observe that this new interpolator generalizes in a nearly identical way to the optimal response-linear achievable interpolator.
Finally, we show that the concept of optimal response-linear achievable interpolation can be extended to more complex models by providing analogous results for a random features model under the same linear data-generating regime as also considered in (Mei and Montanari, 2019), for instance.

heading: 2 Problem setup
text: In this paper we investigate overparametrized linear regression. We assume there exists w? ∈ Rd (unknown) so that yi = 〈w?, xi〉+ ξi for i ∈ {1, . . . , n}, with i.i.d. noise ξi ∈ R (unknown) such that E(ξi) = 0,E(ξ2i ) = σ2 and i.i.d. features xi ∈ Rd that follow a distribution Px with mean E(xi) = 0 and covariance matrix E(xixTi ) = Σ. We store the features in a random matrix X ∈ Rn×d with rows xi ∈ Rd, the response variable in a random vector y ∈ Rn with entries yi ∈ R, and the noise in a random vector ξ ∈ Rn with entries ξi ∈ R. Throughout the paper we assume that d ≥ n. We consider the whitened data matrix Z = XΣ−
1 2 ∈ Rn×d, whose rows satisfy E(zizTi ) = Id, where
Id ∈ Rd×d is the identity matrix. We place a prior on the true parameter in the form w? ∼ Pw? such that E(w?) = 0 and E(w?w?T ) = r 2
d Φ. Here, Φ is a positive definite matrix and r 2 is the
signal. We sometimes abuse terminology and refer to Φ as the covariance of the prior even though r2
d Φ is the covariance matrix. Our results will be proved in general, but for the sake of exposition it can be assumed that xi ∼ N (0,Σ), ξ ∼ N (0, σ2In) and w?∼N (0, r 2
d Φ). We also define the signal-to-noise ratio δ = r2/σ2 and consider the squared error loss ` : (x, y) ∈ R2 7→ (x − y)2. Througout the paper, we assume the following two technical conditions hold.
Assumption 1. Px(xi ∈ V ) = 0 for any linear subspace V of Rd with dimension smaller than d. Assumption 2. For all Lebesgue measurable sets A ⊆ Rd, ν(A) > 0 implies Pw?(w? ∈ A) > 0, where ν is the standard Lebesgue measure on Rd.
Assumption 1 is needed only so that rank(X) = n with probability 1 (for a proof see A.4). A sufficient condition is that Px has a density on Rd. A sufficient condition for Assumption 2 is that Pw? has a positive density on Rd. Now, our goal is to minimize the population risk
r(w) = Ex̃,ξ̃ ( (〈w, x̃〉 − ỹ)2 ) ,
or, equivalently, the excess risk r(w)− r(w?). Here, (x̃, ỹ, ξ̃) is a random variable which follows the distribution of (x1, y1, ξ1), . . . , (xn, yn, ξn) and is independent from them. Throughout the paper we write Ezg(z, z̃) to denote the conditional expectation E(g(z, z̃)|z̃), for two random variables z and z̃ and for a function g. The population risk satisfies
r(w) = (w − w?)TΣ(w − w?) + σ2 = ‖w − w?‖2Σ + r(w?), (1)
where ‖w‖2Σ = wTΣw. We define the bias and variance of w ∈ Rd by the decomposition
Eξ,w? r(w) = B(w) + V (w), (2)
where
B(w) = Eξ,w?‖E(w|w?, X)− w?‖2Σ V (w) = Eξ,w?‖w − E(w|w?, X)‖2Σ. (3)
One of the main paradigms to minimize the (unknown) population risk is based on minimizing the empirical risk R(w) = 1n ∑n i=1(〈w, xi〉 − yi)2 = 1 n ∑n i=1 `(w
Txi, yi) (Vapnik, 1995). In our setting, minimizing the empirical risk is equivalent to finding w ∈ Rd such that Xw = y.

heading: 3 Interpolators
text: An interpolator is any minimizer of the empirical risk. Let G be the set of interpolators, which in linear regression can be written as
G = {w ∈ Rd : Xw = y}.
As rank(X) = n with probability 1, we have G 6= ∅ with probability 1. In linear regression, the implicit bias of gradient descent initialized at 0 is the minimum-norm interpolator (Gunasekar et al., 2018). We define the minimum-norm interpolator by
w`2 = arg min w∈Rd :Xw=y
‖w‖22 = X†y,
where X† ∈ Rn×d is the Moore-Penrose pseudoinverse (Penrose, 1955). As rank(X) = n, we can also write X† = XT (XXT )−1. The second interpolator of interest is a purely theoretical device, previously used in (Muthukumar et al., 2019) to specify a fundamental limit to how well any interpolator in linear regression can generalize.
Definition 1. The best possible interpolator is defined as
Wb = arg min w∈G r(w).
We can write
Wb = arg min w∈Rd :Xw=y
‖Σ 12 (w − w?)‖22,
and after a linear transformation and an application of a result on approximate solutions to linear equations (Penrose, 1956), we obtain
Wb = w ? + Σ− 1 2 (XΣ− 1 2 )†ξ. (4)
We notice that the best possible interpolator fits the signal perfectly by having access to the true parameter w? and fits the noise through the term Σ− 1 2 (XΣ− 1 2 )†ξ by having access to the noise vector ξ in the training data. In general, this interpolator cannot be implemented as it requires access to the unknown quantities w? and ξ. We are interested in interpolators which can be achieved by some algorithm using the data X and y.
Definition 2. We define an estimator w ∈ Rd to be achievable if there exists a function f such that w = f(X, y,Σ,Φ, δ).
In our definition of achievability, we allow for knowledge of the population data covariance, the signal-to-noise ratio, and the prior covariance to define a fundamental limit to what generalization performance can be achieved also without access to these quantities, and we later empirically show that we can successfully approach this limit using only the knowledge of the training data X and y, in considered examples (see Section 5.1). Moreover, our theory is also useful in situations when one has access to some prior information about the regression problem which they can incorporate into an estimate of Σ, δ, Φ (for example, one may know the components of xi are independent and hence Σ is diagonal) and hence it is relevant to consider a broader class than w = f(X, y).
Definition 3. We define the set of response-linear achievable estimators by
L = { w ∈ Rd : ∃f such that w = f(X, y,Σ,Φ, δ) where y ∈ Rn 7→ f(X, y,Σ,Φ, δ) is linear }.
Linearity of y ∈ Rn 7→ f(X, y,Σ,Φ, δ) is equivalent to f(X, y,Σ,Φ, δ) = g(X,Σ,Φ, δ)y, where g is any function which has image in Rd×n. The notion of optimality that we introduce is that of the optimal response-linear achievable interpolator, which is the interpolator that minimizes the expected risk in the class L.
Definition 4. We define the optimal response-linear achievable interpolator by
wO = arg min w∈G∩L
Eξ,w?r(w)− r(w?). (5)

heading: 4 Main results
text: By definition, the interpolator wO has the smallest expected risk among all response-linear achievable interpolators. Our first contribution is the calculation of its exact form.
Proposition 1. The optimal response-linear achievable interpolator satisfies
wO =
( δ
d ΦXT +Σ− 1 2 (XΣ−
1 2 )† )( In+ δ
d XΦXT
)−1 y. (6)
For an isotropic prior Φ = Id, wO depends only on the population covariance Σ and the signal-tonoise ratio δ so that wO can be approximated using estimators of these quantities, which is what we do in Sections 5 and A.8. Even if Φ 6= Id, one might have some information about the prior covariance, which can be incorporated into an estimate Φ̂ and used instead of Φ. However, even if no such estimate is available, in Section A.9 we empirically show that, in our examples, using Φ̂ = Id when Φ 6= Id has a small effect on generalization. Secondly, using results of Gunasekar et al. (2018) on the implicit bias of converging mirror descent, we show that the optimal response-linear interpolator is the limit of gradient descent preconditioned by the inverse of the population covariance, provided that it converges and is suitably initialized.
Proposition 2. The optimal response-linear achievable interpolator is the limit of preconditioned gradient descent wt+1 = wt − ηtΣ−1∇R(wt), (7) provided that the algorithm converges, initialized at
w0 = δ
d ΦXT
( In + δ
d XΦXT
)−1 y. (8)
The interpolator wO does not have the smallest bias or the smallest variance in the bias-variance decomposition Eξ,w? r(w) = B(w) + V (w), but rather achieves a balance. This is related to the results of (Amari et al., 2021). Their setting looks at interpolators achieved as the limit of preconditioned gradient descent in linear regression (preconditioned with some matrix P ) and initialized at 0. Such interpolators can be written as w = PXT (XPXT )−1y. For these interpolators, they compute the risk of w, separate the risk into a variance and a bias term and using random matrix theory they find what the variance and bias terms converge to when d → ∞, n → ∞ in a way such that d/n → γ > 1. For these calculations to hold, they assume that the spectral distribution of (Σd)d∈N converges weakly to a distribution supported on [c, C] for some c, C > 0. Then, after obtaining the limiting variance and bias, they prove which matrices P minimize these limits separately (not their sum, which is the overall asymptotic risk).
We approach the problem from the other direction. That is, we do not a priori consider interpolators that can be achieved as limits of specific algorithms, but we directly look at which interpolator minimizes the risk as a whole (not bias and variance separately). Only after computing the optimal response-linear interpolator, we show in Proposition 2 that the interpolator is in fact the limit of preconditioned gradient descent, however with a specific initialization. Our results hold for every finite d ≥ n and we do not put assumptions on the eigenvalues or the spectral distribution of Σ. In particular, we can recover the results of (Amari et al., 2021) as a special case of Proposition 1. If we take the signal-to-noise ratio δ → 0 (by taking r2 → 0) in Proposition 1, we obtain the matrix P
which achieves optimal variance and if we take δ →∞ (by taking σ2 → 0), we obtain the matrix P which achieves optimal bias. Moreover, we provide a further extension in Proposition 3.
We show that the preconditioned gradient descent wt+1 = wt − ηtΣ−1∇R(wt) achieves optimal variance among all interpolators when initialized at any deterministic w0 and for any finite d, n ∈ N.
Proposition 3. The limit of preconditioned gradient descent wt+1 = wt−ηtΣ−1∇R(wt) initialized at a deterministic w0 ∈ Rd, provided that it converges, satisfies
lim t→∞ wt = arg min w∈G V (w). (9)
We note that the optimal variance is achieved among all interpolators, not only among response-linear achievable interpolators.
A natural question to ask is whether the optimal response-linear achievable interpolator wO provides a significant benefit compared to other interpolators. A second question is whether we can successfully approximate the optimal response-linear achievable interpolator without knowledge of the population covariance Σ and the signal-to-noise ratio δ. In the following section, we illustrate that both the interpolator with optimal variance and the interpolator with optimal bias can generalize arbitrarily badly in comparison to wO as a function of the eigenvalues of the population covariance. In the same regimes where this happens, we present numerical evidence that we can successfully approximate wO by an empirical interpolator wOe without any prior knowledge of Σ or δ by using the Graphical Lasso estimator (Friedman et al., 2007) of the covariance matrix Σ and choosing the empirical estimate of δ by crossvalidation on a subset of the data.

heading: 5 Comparison of interpolators
text: First, we present an example where the minimum-norm interpolator w`2 generalizes arbitrarily worse than the best response-linear achievable interpolator wO. Second, we give an example where an interpolator with optimal variance generalizes arbitrarily worse than wO. This shows that arbitrarily large differences in test error are possible within the class of estimators with zero training error.
In the examples, we consider a setting where xi ∼ N (0,Σ) and w? ∼ N (0, r 2
d Φ). Therefore, throughout Section 5 we assume Px = N (0,Σ) and Pw? = N (0, r 2
d Φ). Before presenting these examples, we discuss approximating wO by an interpolator, wOe, which uses only the data X and y.

heading: 5.1 Empirical approximation
text: The interpolator wO is the limit of the algorithm wt+1 = wt − ηtΣ−1∇R(wt). (10)
The population covariance Σ is required to run this algorithm. However, the matrix Σ is usually unknown in practice. One may want to estimate Σ. However, if one replaces Σ by Σ̃ = XTX/n+λId (with λ ≥ 0), then the limit of (10) is the same as the limit of gradient descent (provided that both algorithms converge). This is because, using the singular value decomposition of X , one can show
Σ̃− 1 2 (XΣ̃− 1 2 )†y = X†y.
The preconditioned gradient update wt+1 − wt = ηtP−1∇R(wt) has to not belong to Im(XT ) in order to not converge to the minimum-norm interpolator. Hence, using P = Σ̃ = XTX/n+λId (for example, also the Ledoit-Wolf shrinkage approximation (Ledoit and Wolf, 2004)) in preconditioned gradient descent removes the benefit of preconditioning in terms of generalization of the limit.
We use the Graphical Lasso approximation (Friedman et al., 2007). We empirically observe that in the examples considered in this paper (Figures 1, 2, 3, 4, 5, 6) using the Graphical Lasso covariance Σe instead of Σ has nearly no effect on generalization. Under specific assumptions, Ravikumar et al. (2011) provide some convergence guarantees of the Graphical Lasso.
In regards to approximating the signal-to-noise ratio δ, we choose δe that minimizes the crossvalidated error on random subsets of the data. In this way, we arrive at the interpolator
wOe= ( δe d XT +Σe − 12 (XΣe − 12 )† )( In+ δe d XXT )−1 y, (11)
which approximates wO and is a function of only X and y. We note that the interpolator wOe uses Id in place of the prior covariance matrix.
In the experiments (Figures 1, 2, 3, 4, 5, 6) we used the Graphical Lasso implementation of scikitlearn (Pedregosa et al., 2011) with parameter α = 0.25 (α can also be crossvalidated for even better performance) and in estimating δ, for each δe in {0.1, 0.2, . . . , 1, 2, . . . , 10}, we computed the validation error on a random, unseen tenth of the data and averaged over 10 times. The δe with smallest crossvalidated error was chosen.

heading: 5.2 Random matrix theory concepts
text: For presenting the discussed examples we need to review some concepts from random matrix theory.
Definition 5. For a symmetric matrix Σ ∈ Rd×d with eigenvalues λ1 ≥ λ2 ≥ · · · ≥ λd ≥ 0 we define its spectral distribution by FΣ(x) = 1d ∑d i=1 1[λi,∞)(x).
The following assumptions will be occasionally considered for the covariance matrix Σ.
Assumption 3. There exists kmax > 0 such that λmax(Σ) ≤ kmax uniformly for d ∈ N. Assumption 4. There exists kmin > 0 such that kmin ≤ λmin(Σ) uniformly for d ∈ N. Assumption 5. The spectral distribution FΣ of the covariance matrix Σ converges weakly to a distributionH supported on [0,∞).
Marčenko and Pastur (1967) showed that there exists a distribution F̃γ such that FZT Z n −→ F̃γ , weakly, with probability 1 as n→∞, d→∞ with d/n→ γ. In our discussion, xi = Σ 1 2 zi, where zi∼N (0, Id) independently. Then, under Assumption 5, it can be shown that the spectral distribution of Σ̂ = XTX/n = Σ 1 2ZTZΣ 1 2 /n converges weakly, with probability 1 to a distribution supported on [0,∞), which we denote by Fγ , see e.g. (Silverstein and Choi, 1995). Similar arguments also show that the spectral distribution of XXT /n ∈ Rn×n converges weakly, with probability 1.
Definition 6. For a distribution F supported on [0,∞), we define the Stieltjes transform of F , for any z ∈ C \ R+ by
mF (z) = ∫ ∞ 0 1 λ− z dF(λ).
The weak convergence of the spectral distribution of Σ̂ to Fγ is equivalent to mΣ̂(z)→ m(z) and mXXT /n(z)→ v(z) almost surely for all z ∈ C \ R+, where m and v are the Stieltjes transforms of Fγ and the limiting spectral distribution of XXT /n, respectively (see e.g. Proposition 2.2 of (Hachem et al., 2007)). We call v the companion Stieltjes transform of Fγ .

heading: 5.3 Diverging variance of interpolator with optimal bias
text: Using w0 = 0 in Proposition 3, we choose the interpolator with optimal variance to be (see A.3)
wV = Σ − 12 (XΣ− 1 2 )†y. (12)
When Φ = Id, the interpolator with best bias among response-linear achievable interpolators is the minimum-norm interpolator (see Section A.5). We identify an example where the minimum-norm interpolator w`2 generalizes arbitrarily worse than the best response-linear achievable interpolator wO. For this, we exploit results of Hastie et al. (2019) on computing the asymptotic risk of the minimum-norm interpolator. They show that if Φ = Id, under Assumptions 3, 4, 5 and if dn → γ > 1 with n→∞, d→∞ then with probability 1,
Eξ,w?r(w`2)− r(w?) −→ r2
γv(0) + σ2
( v′(0) v(0)2 − 1 ) , (13)
where v is the companion Stieltjes transform introduced in Section 5.2. In comparison, similarly as in (Amari et al., 2021), the asymptotic risk of the best variance estimator wV satisfies that under
Assumption 3 and 5, if n, d→∞ with dn → γ > 1 then with probability 1 we have that
lim d→∞ Eξ,w∗r(wV )− r(w?) = r2 γ − 1 γ ∫ ∞ 0 s dH(s) + σ 2 γ − 1 . (14)
An alternative way is to write ∫∞
0 s dH(s) = limd→∞ Tr(Σ). This result follows by an application
of Theorem 1 of (Rubio and Mestre, 2011), which is in the supplementary material for completeness.
Now we find a regime of covariance matrices Σ, for which the variance term of the minimumnorm solution, V`2 = σ 2( v ′(0) v(0)2 − 1), diverges to infinity, while the risk of wO stays bounded and close to optimal. For this, we consider a generalization of the spike model of covariance matrices (Baik and Silverstein, 2006; Johnstone, 2001), which is a fundamental model in statistics. Here Σ = diag(ρ1, . . . , ρ1, ρ2, . . . , ρ2) ∈ Rd×d, where the number of ρ1s is d · ψ1 with ψ1 ∈ [0, 1]. This model was also considered in (Richards et al., 2021) where it is called the strong weak features model. In this regime, it is possible to explicitly calculate the companion Stieltjes transform v(0) and v′(0) of (13). In the case that γ = 2, ψ1 = 1/2 we have
V`2 = σ 2
( v′(0) v(0)2 − 1 ) = σ2 2 (√ ρ1 ρ2 + √ ρ2 ρ1 + 2 ) . (15)
If we fix ρ1 = 1 and take ρ2 → 0, then the variance term V`2 diverges to infinity. This also means that the asymptotic risk of the minimum-norm interpolator diverges to infinity. Moreover, the asymptotic risk of wV in (14) evaluates to
lim d→∞
Eξ,w∗r(wV )− r(w?) = ( ψ1ρ1 + (1− ψ1)ρ2 )( 1− 1
γ
) + σ2
γ − 1 . (16)
In addition, by construction of wO, we know that Eξ,w?r(wO) ≤ Eξ,w?r(wV ) and therefore the asymptotic limit of Eξ,w?r(wO)− r(w?), as d/n→ γ > 1, stays bounded by (16) as ρ2 → 0. The expected generalization error in the setting described above is illustrated in Figure 1.
We note that the empirical estimator wOe (yellow points), which is a function of only the training data X and y and does not use the population covariance Σ or the signal-to-noise ratio δ, performs almost identically to the optimal response-linear achievable interpolator wO (cyan points).
In this example, we chose γ = 2 and ψ1 = 1/2 deliberately. One does not achieve diverging variance for an arbitrary choice of γ and ψ1. However, for any γ > 1 such that γψ1 = 1, the phenomenon of Figure 1 holds (see A.7 of the supplementary material).

heading: 5.4 Diverging bias of interpolator with optimal variance
text: Now, we illustrate a regime where the best variance interpolator wV generalizes arbitrarily worse than wO. In the same strong and weak features covariance model described above in Section 5.3, when γ = 2 and ψ1 = 1/2, if we instead have ρ1 → ∞ and ρ2 = 1, then the asymptotic risk (16) diverges to infinity linearly. However, the variance of the minimum-norm interpolator in (15) diverges only like √ ρ1. Moreover, the bias term satisfies
B`2 = r2 γv(0) = r2 γ √ ρ1ρ2,
which also diverges like √ ρ1. Now, because Eξ,w?r(wO) ≤ Eξ,w?r(w`2), we have that
lim d→∞
Eξ,w?r(wO) ≤ r2
γ
√ ρ1ρ2 + σ2
2 (√ ρ1 ρ2 + √ ρ2 ρ1 + 3 ) ,
so that the asymptotic risk of wO diverges to infinity as √ ρ1. We illustrate this in Figure 2.
We notice that the empirical approximation wOe again performs in a nearly identical way to the optimal response-linear achievable interpolator wO. Moreover, importantly, we note that wV and wO are limits of the same algorithm, wt+1 = wt − ηtΣ−1∇R(wt), only with different initialization. Hence, this shows that different initialization of the same optimization algorithm can have an arbitrarily large influence on generalization through implicit bias.

heading: 6 Random features regression
text: The concept of optimal interpolation as a function which is linear in the response variable, is general and can be extended beyond linear models. We present an extension of Proposition 1 to the setting of random features regression. Random features models were introduced as a random approximation to kernel methods (Rahimi and Recht, 2008) and can be viewed as a two-layer neural network with first layer randomly initialized and fixed as far as training is concerned. They can be shown to approximate neural networks in certain regimes of training and initialization and hence are often considered in the literature as a first step to address neural networks (e.g. (Jacot et al., 2018)). We consider data generated in the same way as before, yi = 〈xi, w?〉 + ξi, and the model to be a two-layer neural network fa : Rd 3 x 7→ aTσ(Θx/ √ d), where the first layer Θ ∈ RN×d is randomly initialized. This
setting, along with xi and rows of Θ belonging to the sphere Sd−1( √ d) with radius √ d in Rd, is often considered in the literature on interpolation of random features models (Mei and Montanari, 2019; Ghorbani et al., 2021). If we analogously define the optimal response-linear achievable interpolator in random features regression by
aO = arg min a∈G∩L
Eξ,w?r(fa)− r(w?), (17)
where here G = {a ∈ RN : Za = y} is the set of interpolators, Z = σ(XΘT / √ d) and L is the same as in Definition 3, then the following analogue of Proposition 1 holds.
Proposition 4. The optimal response-linear achievable interpolator (17) in random features regression satisfies
aO=Σ −1 z ( ΣzxΦX T +ZT ( ZΣ−1z Z T )−1(d
δ In+XΦX
T−ZΣ−1z ΣzxΦXT ))(d
δ In+XΦX
T )−1 y,
Here Σz = Ex̃(σ(Θx̃/ √ d)σ(Θx̃/ √ d)T ) and Σzx = Ex̃(σ(Θx̃/ √ d)x̃T ) are covariance and crosscovariance matrices, respectively. This interpolator can be again obtained as the implicit bias of preconditioned gradient descent using results of Gunasekar et al. (2018).
Proposition 5. The optimal response-linear achievable interpolator (17) in random features regression is the limit of preconditioned gradient descent on the last layer,
wt+1 = wt − ηtΣz−1∇R(wt), provided that the algorithm converges, initialized at
a0 = Σ −1 z ΣzxΦX T
( d
δ In +XΦX
T )−1 y.
In Section A.12, we illustrate the test error of fa, with a = aO in comparison to the test error for the minimum-norm interpolator a = a`2 = Z †y on a standard example.

heading: 7 Conclusion
text: In this paper, we investigated how to design interpolators in linear regression which have optimal generalization performance. We designed an interpolator which has optimal risk among interpolators that are a function of the training data, population covariance, signal-to-noise ratio and prior covariance, but does not depend on the true parameter or the noise, where this function is linear in the response variable. We showed that this interpolator is the implicit bias of a covariance-based preconditioned gradient descent algorithm. We identified regimes where other interpolators of interest are arbitrarily worse using computations of their asymptotic risk as dn → γ > 1 with d, n→∞. In particular, we found a regime where the variance term of the minimum-norm interpolator is arbitrarily large compared to our interpolator. This confirms the phenomenon that implicit bias has an important influence on generalization through the choice of optimization algorithm.
We identified a second regime where the interpolator that has best variance is arbitrarily worse than our interpolator. In this second example, both interpolators are the implicit bias of the same algorithm, but with different initialization. This contributes to illustrating that initialization has an important influence on generalization.
We also considered an empirical approximation of the optimal response-linear achievable interpolator, which uses only the training data X and y and does not assume knowledge of the population covariance matrix, the signal-to-noise ratio or the prior covariance and empirically observe that it generalizes in a nearly identical way to the optimal response-linear achievable interpolator in the examples that we consider.
A limitation of this work includes a precise guarantee on the approximation error of the Graphical Lasso for a general covariance matrix Σ. Some guarantees are in (Ravikumar et al., 2011), however establishing guarantees for a general covariance matrix would be a contribution on its own.
A natural question for future research, which also motivated our work, is how to systematically design new ways of interpolation, which are adapted to the distribution of the data and related to notions of optimality, for more complex overparametrized machine learning models such as neural networks.

heading: 8 Acknowledgements
text: The authors would like to thank Dominic Richards, Edgar Dobriban and the anonymous reviewers for valuable insights which contributed to the technical quality of the paper. Eduard Oravkin was part-time employed at the Department of Statistics during a part of this project. Patrick Rebeschini was supported in part by the Alan Turing Institute under the EPSRC grant EP/N510129/1.

Review summaries:
This paper studies the generalization error of a class of interpolator in linear regression and identifies conditions for optimality. They consider a setting where there is a true linear model generating the data, whose parameters are drawn from a Gaussian distribution with zero-mean and covariance Phi. In this setting they consider a natural class of estimators that are a function of the data, the population covariance matrix, the prior covariance matrix Phi and the signal to noise ratio. Among this class of estimators they calculate a closed-form expression for the optimal interpolator, with respect to the expected excess risk and show that there are certain settings where the expected excess risk of this estimator can be much lower than the risk of the OLS solution. They also provide an extension of their results to the random features model.

The paper under reviews studies the problem of optimal linear interpolation. More precisely the authors derive an explicit formula for the *best* linear interpolator that can be written as a linear function of the outputs and as a function of desired quantities. From this, they show how to recover theoretically this estimator from a optimization algorithm. Finally they expose cases where the designed estimator show some good generalization properties while the miminum norm interpolator does not.  The link with kernel methods is also made.

This paper is concerned with inducing a certain type of best possible interpolator in linear regression. It is based on the notion of best possible interpolator, defined in Muthukumar et al., 2019., which is the best estimator (i.e., with minimum risk) among the class of interpolators. Based on that, the authors define a new theoretical estimator, the best estimator among interpolators that are also linear in the response variable.  They manage to induce a closed form of this estimator and prove that it is the limit of preconditional gradient descent under a specific initialization. In addition, using approximations for quantities that are usually unknown in practice (e.g.  population covariance matrix, signal to noise ratio), they turn this theoretical estimator into a fully empirical one and compare its performance to other frequently used estimators, showing that this estimator performs arbitrarily better in certain cases. Finally, they extend this notion of best possible theoretical estimator to the random features setting, where they again induce a closed form for this estimator and prove that it is the limit of preconditioned gradient descent under a specific initialization. 

This paper considers the problem of optimal interpolation in linear regression. It defines the class of response-linear achievable interpolators, and provides a closed-form expression for the optimal interpolator in this class. The paper is about an important topic in machine learning and statistics. However, it is unclear to me whether the "achievable" interpolators are really achievable, and as a result if the investigation is meaningful.