# Flat Minima in Linear Estimation and an Extended Gauss Markov Theorem

Simon Segert
Princeton Neuroscience Institute
Princeton University
Princeton, NJ
<EMAIL_ADDRESS>

###### Abstract

We consider the problem of linear estimation, and establish an extension of
the Gauss-Markov theorem, in which the bias operator is allowed to be non-zero
but bounded with respect to a matrix norm of Schatten type. We derive simple
and explicit formulas for the optimal estimator in the cases of Nuclear and
Spectral norms (with the Frobenius case recovering ridge regression).
Additionally, we analytically derive the generalization error in multiple
random matrix ensembles, and compare with Ridge regression. Finally, we
conduct an extensive simulation study, in which we show that the cross-
validated Nuclear and Spectral regressors can outperform Ridge in several
circumstances.

##  1 Introduction

Linear models are among the most used of all machine-learning models in
applications to science and engineering. In addition to this practical
interest, they are also of great theoretical interest. Indeed, the empirical
successes of neural networks have proven somewhat at odds with traditional
received wisdom from classical statistical learning theory(Belkin et al.,
2019). This has begun to be reconciled by careful analysis of well-chosen
linear “model systems” such as high-dimensional regression (Hastie et al.,
2019), kernel ridge regression (Canatar et al., 2021; Jacot et al., 2020), or
linear neural networks (Saxe et al., 2013), which can reproduce many
qualitative properties of learning dynamics of neural networks. Thus, the
careful study of linear models, especially in the limit of a large number of
features and observations can prove highly valuable for qualitative
understanding of non-linear models.

Additionally, it has been found that explicit regularization is often a key
ingredient for achieving high performance of both linear and non-linear
models. For example, neural networks are often trained using L2 regularization
(Goodfellow et al., 2016) or dropout (Srivastava et al., 2014). By contrast,
most of the above-mentioned theoretical work on linear models considers
specifically L2 regularization (i.e. Ridge regression). It is quite natural,
then, to think that the linear setting could be used as a test case for
studying other forms of regularization. This question becomes especially
interesting in light of several recent works which have explored the effect of
various non-standard forms of regularization on neural networks such as rank
constraints (Kamalakara et al., 2022; Yang et al., 2020) or Spectral norms
(Johansson et al., 2022; Yoshida & Miyato, 2017). The Nuclear norm (or
convexified rank) also plays a key role in other kinds of high-dimensional
non-linear learning problems (Hu et al., 2021), most notably matrix
completion, and is related to dropout regularization in neural networks
(Mianjy & Arora, 2019).

Thus with the general motivation of gaining a detailed understanding of
different kinds of regularizations in a tractable yet previously-validated and
informative setting, we study in this work the effect of different matrix norm
regularizations in the context of high-dimensional linear regression. More
specifically, we consider the problem of linear regression, under the
constraint that the bias matrix is bounded according to some fixed matrix
norm. We prove a Gauss-Markov-like theorem for this setting, which
characterizes the optimal estimator at a fixed level of bias, assuming the
norm in question is of Schatten type. We will be specifically interested in
Nuclear, Frobenius, and Spectral norms, in which case the optimal estimators
have especially simple forms, which we derive.

Next, we characterize the test error of Nuclear and Spectral estimators in
several random matrix ensembles (the corresponding results for Ridge
regression being well known). We find an intriguing pattern: while the Ridge
estimator can always attain the lowest error of any estimator, the minima can
be sharp (as a function of regularization strength). The Nuclear estimator, by
contrast, can often attain test error nearly as low as the Ridge case, but
with significantly flatter minima. This has implications for practical
situations when the regularization strength is selected using a noisy search
procedure such as cross validation, as it could happen that such a procedure
finds a better solution in the latter case, even if the actual global optimum
is better in the former case. We then perform a series of simulation studies
with cross-validation, in which we find that this can actually happen, with
both Gaussian features and Random Fourier features, and more generally
characterize the relative performance of the three estimators across a range
of several generative factors. We close with a survey of related work and
general discussion.

##  2 Theory

###  2.1 Setup and Extended Gauss-Markov Theorem for Schatten norms

Throughout, we will be concerned with the classic regression problem. Denote
the training data by $X$, which is an $N\times d$ matrix, and the training
targets by $Y$, which is an $N\times 1$ vector. $Y$ is assumed to have
distribution $Y=X\beta_{0}+\epsilon$, where $\epsilon$ are independent samples
from a fixed noise distribution.

We will assume that the noise distribution has mean zero and variance
$\sigma^{2}=1$ , but otherwise do not place distributional assumptions on it.
We will restrict attention to the class of Linear Estimators. These are models
that take the form

$\hat{\beta}=LY$ (1)

Where $L$ is an $d\times N$ matrix that depends (possibly non-linearly) on the
training data. Given a new datapoint $x_{test}$,the predicted value is just
$\hat{y}=\langle x_{test},\hat{\beta}\rangle$. One simple observation is that
for any linear estimator the bias $\mathbb{E}_{\epsilon}\hat{\beta}-\beta_{0}$
can be expressed as a linear function of $\beta_{0}$:

$\mathbb{E}_{\epsilon}\hat{\beta}-\beta_{0}=(LX-I)\beta_{0}$ (2)

We thus define the bias operator $B:=LX-I$. Note that the actual bias value
generally depends on the unknown quantity $\beta_{0}$, while the variance
$Var_{\epsilon}(\hat{\beta})=LL^{T}$ does not. The classic Gauss-Markov
theorem tells us that if we want to impose exactly zero bias111Note the subtle
point that this is the only case in which we can exactly control the bias
without knowing the value of $\beta_{0}$, then the minimal variance estimator
is the ordinary least squares estimator. But what if we want to relax this, to
allow non-zero but bounded amount of bias? Taking this desideratum literally,
we run into the fundamental issue that the bias depends on the unknown
$\beta_{0}$. As a proxy, we propose to control the size of $B$ instead, using
some choice of matrix norm. In what follows, we will only consider Schatten
norms, $\|M\|_{p}=(\sum_{i}\sigma_{i}(M)^{p})^{1/p}$, $p\geq 1$, although in
principle other matrix norms could be used. 222We reserve the notation
$\|\cdot\|_{p}$ for a Schatten norm of a matrix, and use $|\cdot|_{p}$ for
Euclidean vector normsThis motivates the following:

######  Definition 1.

Let $p\geq 1$, and let $C\geq 0$. For any matrix $X\in\mathbb{R}^{N\times
d}$,the p-Bias constrained Linear estimator $L_{p}(X)$ is defined as the
minimal variance estimator which has p-bias of at most $C$. That is,

$L_{p}(X)=argmin_{L\in\mathbb{R}^{d\times N};\|LX-I_{d}\|_{p}\leq
C}Tr(LL^{T})/2$ (3)

where $\|\cdot\|_{p}$ is a Schatten norm. If $Y$ is a vector of regression
targets, the estimated coefficient vector is $\hat{\beta}:=L_{p}(X)Y$.

We now state:

######  Theorem 2.

(Extended Gauss-Markov) Let $X\in\mathbb{R}^{N\times d}$ be a matrix of
observations. Assume that $G:=X^{T}X$ is invertible, and let
$G=Udiag(\sigma)U^{T}$ be the diagonalized form. For $p\geq 1$, the p-Bias
constrained estimator with bound $C$ can be expressed in the form
$L_{p}(X)=\hat{G}^{-1}X^{T}$, where $\hat{G}$ is symmetric,simultaneously
diagonalizable with $G$, and satisfies $G\preceq\hat{G}$333i.e., $\hat{G}-G$
is non-negative definite. For the following special cases of $p$ we further
have:

$p=1$ (Nuclear norm): $\hat{G}=Udiag(\max(\sigma,\alpha))U^{T}$

$p=2$ (Frobenius norm): $\hat{G}=G+\alpha I_{d}$

$p=\infty$ (Spectral norm): $\hat{G}=(1+\alpha)G$

where $\alpha\geq 0$ is determined by $C$, and $max$ denotes the elementwise
maximum, i.e. $max(\sigma,\alpha)_{i}=max(\sigma_{i},\alpha)$

Note that the $p=2$ case is just Ridge regression, while the $p=\infty$ case
is just scalar shrinkage of $\hat{\beta}_{OLS}$ towards 0, as per Stein
(1962). We refer to the $p=1$ case as Nuclear Regression and likewise Spectral
Regression for the $p=\infty$ case. The exact relation between $C$ and
$\alpha$ is given in A.8.

We also note that our theorem does not perfectly generalize the classical
Gauss-Markov theorem. Indeed, the classical theorem is usually stated in the
form $Var(\hat{\beta}_{OLS})\preceq Var(\hat{\beta})$ where $\hat{\beta}$ is
any unbiased linear estimator (Wooldride, 2015). However, taking $C=0$ in our
Theorem 2, we would conclude only that $Tr(Var(\hat{\beta}_{OLS}))\leq
Tr(Var(\hat{\beta}))$, which is strictly weaker. But, we regard this as a bit
of a technicality, and think our theorem faithfully captures the spirit of the
original.

Figure 1: Illustration of Theorem 2 Each point on the curve corresponds to a
regression fit with a certain value of $\alpha$, on the matrix
$X=Diag(\sqrt{1},\sqrt{2},\dots,\sqrt{10})$. We show the norm of the bias
matrix $B=LX-I$ on the x axis, and the sample variance $Tr(L^{T}L)$ on the y
axis. The green curve always lies below the other two on the left plot, and
analogously for the other two plots.

###  2.2 Formulas for Test Error

We next analytically characterize the generalization performance of the
estimators derived in the previous section, under certain specific assumptions
on the distribution of $X$. We will consider two different random matrix
ensembles: a standard spherical Gaussian ensemble, and another with orthogonal
predictors but arbitrary spectrum. We will compute the test error averaged
over random datasets, in an appropriate thermodynamic limit.

####  2.2.1 Spherical Gaussian ensemble

Consider a predictor matrix $X\in\mathbb{R}^{N\times d}$, where each entry is
an independent centered Gaussian of variance $1/N$. We consider the
thermodynamic limit in which $N,d\to\infty$ and $d/N\to\lambda<1$. The
training targets are distributed as $Y_{i}=\langle
X_{i},\beta_{0}\rangle+\sigma\epsilon_{i}$, where $\epsilon_{i}\sim N(0,1)$,
and the magnitude of the ground-truth coefficient vector $\beta_{0}$ satisfies
$|\beta_{0}|^{2}/d\to\beta^{2}$, for some fixed $\beta>0$.

The testing error for a given dataset $X,Y$ is defined as

$MSE:=\mathbb{E}_{x\sim N(0,I_{d}/N)}(\langle x,\beta_{0}\rangle-\langle
x,\hat{\beta}\rangle)^{2},$ (4)

where $\hat{\beta}$ is the coefficient vector estimated from $X$ and $Y$.
(Note that we do not include exogenous noise in the testing data; if we did,
it would just contribute an additive offset of $\sigma^{2}$). We observe that
the expression can be simplified to $|\beta_{0}-\hat{\beta}|^{2}/N$, however
we write it in the more verbose form above to make clear that it is obtained
as an average over the distribution of the testing point $x$.

We are interested in the average value of this over different datasets. Using
standard techniques in Random Matrix theory, it is not hard to show (cf. A.6):

######  Proposition 2.1.

Consider the $p$-bias constrained estimator with regularization strength
$\alpha$.

In the above limit, the average test error
$Err_{p}(\alpha):=\mathbb{E}_{X,Y}MSE$ is given by

$Err(\alpha)=\lambda\int\beta^{2}(1-{\frac{x}{f_{\alpha}(x)}})^{2}+\sigma^{2}{\frac{x}{f_{\alpha}(x)^{2}}}\mu^{\lambda}_{MP}(dx)$
(5)

where $\mu_{MP}^{\lambda}$ is the Marchenko-Pastur density with concentration
parameter $\lambda$, and $f_{\alpha}$ is defined casewise as follows:

$f_{\alpha}(x)=max(x,\alpha),p=1$ (6) $f_{\alpha}(x)=x+\alpha,p=2$ (7)
$f_{\alpha}(x)=x(1+\alpha),p=\infty$ (8)

In our later analyses, we will employ this formula, and evaluate the integral
using numerical quadrature. However, it is interesting to note that the
integral can be actually be further simplified in each of the three cases. For
the ridge case, there is a well-known formula in terms of the Stieltjes
transform of $\mu_{MP}^{\lambda}$ that easily follows from the above
integral(Bai & Silverstein, 2010). In the spectral case, we can derive a
formula using simple calculus:

######  Proposition 2.2.

For the Spectral estimator in the spherical Gaussian setup, the test error in
the thermodynamic limit equals

$Err_{\infty}(\alpha)={\frac{\lambda\beta^{2}\alpha^{2}+{\frac{\lambda}{1-\lambda}}\sigma^{2}}{(1+\alpha)^{2}}}$
(9)

To state the result for the Nuclear case,it is necessary to introduce the
Appel hypergeometric function $F_{1}$(Appell, 1925), which generalizes the
classical Gauss hypergeometric function to two variables.

######  Proposition 2.3.

For the Nuclear estimator in the spherical Gaussian setup, the test error in
the thermodynamic limit can be expressed in terms of the Appel hypergeometric
function $F_{1}$:

$Err_{1}(\alpha)=\begin{cases}{\frac{\sigma^{2}\lambda}{1-\lambda}}&\alpha\leq\lambda_{-}\\\
{\frac{\sigma^{2}\lambda}{1-\lambda}}+\lambda\beta^{2}CDF_{\lambda}(\alpha)+\sqrt{\alpha-\lambda_{-}}\sum_{r\in\\{-1,1,2\\}}c_{r}(\alpha)F_{1}({\frac{3}{2}},1-r,-{\frac{1}{2}},{\frac{5}{2}},1-{\frac{\alpha}{\lambda_{-}}},{\frac{\alpha-\lambda_{-}}{\lambda_{+}-\lambda_{-}}})&\lambda_{-}\leq\alpha\leq\lambda_{+}\\\
{\frac{\lambda(1+\lambda)}{\alpha^{2}}}+{\frac{\lambda\sigma^{2}-2\lambda\alpha\beta^{2}}{\alpha^{2}}}+\lambda\beta^{2}&\alpha\geq\lambda_{+}\end{cases}$

where $\lambda_{\pm}$ are the limits of the Marchenko-Pastur support:
$\lambda_{\pm}=(1\pm\sqrt{\lambda})^{2}$, $CDF_{\lambda}$ is the cumulative
distribution function of $\mu_{MP}$ and $c_{r}(\alpha)$ are certain rational
functions of $\alpha$.

See A.5 for proof, and precise definition of $F_{1}$.

####  2.2.2 Diagonal ensemble

While the spherical Gaussian case is an instructive starting point, the
assumptions are clearly somewhat limiting. In particular, it implies a very
specific functional form of the spectral density, which may not be a good
match to real data. Thus, we analyze another random matrix ensemble which
allows us the flexibility to specify an essentially arbitrary spectral
density. On the one hand, this new ensemble has its own limiting assumptions
that the Gram matrix $X^{T}X$ is almost surely diagonal, and the observations
(i.e. rows of $X$) are no longer independent. But on the other, having another
random matrix model allows us to assess whether qualitative observations about
the spherical Gaussian case generalize to other settings.

We now describe the random matrix ensemble. We first specify some distribution
$\nu$ supported on $[0,1]$, which will be the spectral density. We also
specify some distribution $\nu_{noise}$ supported on a bounded interval of
$\mathbb{R}_{+}$, which will act as multiplicative random noise on the
spectrum; we require that $E_{x\sim\nu_{noise}}x=1$. To generate training and
testing data, we first sample $\lambda_{i}\sim\nu,i=1,\dots,d$ independently,
and $s_{i}\sim\nu_{noise}$. We then set

$X_{tr}=X_{1}diag(\sqrt{\lambda_{i}s_{i}}),X_{test}=X_{2}diag(\sqrt{\lambda_{i}})$
(10)

where $X_{1}$ and $X_{2}$ are drawn independently from the Stiefel manifold
$O(n,d)$ of $N\times d$ orthogonal matrices. That is $X_{i}^{T}X_{i}=I_{d}$.
We observe that, as noted above: 1.) the Gram matrix $G=X_{tr}^{T}X_{tr}$ is
diagonal almost surely, and 2.) the limiting eigenvalue density of $G$ is
given by $\nu$.

The regression targets are then constructed as before:
$Y_{tr}=X_{tr}\beta_{0}+\epsilon,Y_{test}=X_{test}\beta_{0}$. The test error
is defined as

$MSE:=N^{-1}|X_{test}\hat{\beta}-X_{test}\beta_{0}|^{2}$ (11)

Note that this is slightly different from the expression for MSE in the
spherical Gaussian case (Equation 4). The basic reason is that the
observations (i.e. rows) are here not independent. Thus, in contrast to the
spherical Gaussian case, where we could sample one test observation at a time
and average, here we sample a batch of N (dependent) test observations at a
time, and compute the average error over those N.

As before, we will be interested in the value of this error, marginalized over
all of the randomness, and taken in the thermodynamic limit.

######  Proposition 2.4.

In the thermodynamic limit of the above random matrix ensemble, the average
test error $Err_{p}(\alpha)=\mathbb{E}_{X_{tr},X_{test},\epsilon,s}MSE$ can be
expressed as

$Err_{p}(\alpha)=\lambda\int_{0}^{1}\beta^{2}x(1-{\frac{x}{f_{\alpha}(x)}})^{2}+\sigma^{2}{\frac{x^{2}}{f_{\alpha}(x)^{2}}}\nu(dx)$
(12)

where $f_{\alpha}$ has the same form as in Proposition 2.1.

See A.6 for proof. Note the similarity to the result of Proposition 2.1. The
basic reason for this is that both ensembles have a similar rotational
invariance, which allows the MSE to be expressed as a sum over the spectrum of
$X_{tr}$. Besides the integrating measure, the only difference is an extra
factor of $x$ in the Diagonal case.

In our experiments, we specialize to power law distributions
$d\nu(x)/dx\propto x^{\gamma-1}\textbf{1}_{x\in[0,1]}$, as many real
covariance matrices have approximately this structure (Stringer et al., 2019;
Qin & Colwell, 2018; Liu et al., 1997; Ruderman & Bialek, 1994). In this case,
the expressions for the test error can again be simplified to analytic
expressions that do not involve integrals; the details are straightforward and
left to the interested reader444The Ridge case involves the Gauss
hypergeometric function, the other two only require elementary calculus.

####  2.2.3 On the optimal Regularization form

It is not hard to see that the derivations in the above two sections go
through for any estimator of the form $\hat{\beta}=\hat{G}^{-1}X^{T}Y$, where
the eigenvalues of $\hat{G}$ are related to those of $G$ through some fixed
point-wise transform $\lambda_{i}(\hat{G})=f(\lambda_{i}(G))$, and $\hat{G}$
is simultaneously diagonalizable with $G$. It is thus natural to ask what
would be optimal estimator in this class, if we are allowed to select f
arbitrarily? By formally differentiating the integrand with respect to $f$ and
setting equal to zero555this can be justified using results from Calculus of
Variations, it is easy to see that in both of the above random matrix models,
the optimal form of $f$ corresponds to Ridge regression with regularization
strength $\sigma^{2}/\beta^{2}$.

Does this mean that Ridge regression will always perform better in practice
than any other estimator? Not necessarily. For the above estimator is an
oracle estimator in the sense that it requires knowledge of the ground-truth
$\beta$ and $\sigma$, which are typically not known. In practice, we would
estimate the optimal ridge parameter using,e.g., cross-validation, and hope
that our estimate is close to the oracle-optimal value $\sigma^{2}/\beta^{2}$.
In the finite-sample regime, there will be some noise in the cross-validated
errors, so this essentially amounts to trying to minimize $Err(\alpha)$ given
only access to a small number of noisy estimates of the underlying function.
Additionally, there will be “resolution noise” caused by the fact that we will
typically only try to estimate the test error for a small number of values of
$\alpha$. A very rough rule of thumb is that if we sample $n$ values
$\alpha_{i}$, all within some small distance $\delta$ to the minimum, then
$min_{i}Err(\alpha_{i})\approx\mu+{\frac{(\kappa\delta)^{2}}{(n+1)(n+2)}}$,
where $\mu=\min_{\alpha}Err(\alpha)$ is the true minimum, and $\kappa^{2}$ is
the Hessian of $Err(\alpha)$ at the minimum (see A.3 for justification). We
conclude that the “effective minimum” that we would find using cross
validation is penalized by higher values of curvature, especially for small
$n$.

####  2.2.4 Experimental Validation

In order to validate the above formulas, we generated synthetic datasets
according to the random matrix ensembles from sections 2.2.1 and 2.2.2 , and
compared the average test error with the theoretically predicted values. For
each such dataset, we used $N=100$ training observations, and obtained the
empirical test error by averaging over $100$ datasets. We did this for both
the spherical Gaussian case, as well as the diagonal case with a power law
spectral density $d\nu/dx\propto x^{\gamma-1}\textbf{1}_{x\in[0,1]}$. In
Figure 2 we show several plots with different values of the other
hyperparameters; in all cases we see a very close correspondence between the
empirical and predicted errors. In accordance with the discussion in the
previous section, we can also see visually that the minima of the Nuclear
estimator are often noticeably flatter than those of the Ridge, while also
being only moderately higher.

Figure 2: Predicted vs. empirical test errors, for the spherical Gaussian case
(left) and diagonal case with power law spectrum (right). Each dot represents
an average over 20 different datasets, each with $N=100$ training
observations.Note that the minima for the Ridge curve are notably steeper than
for the other two.

##  3 Experiments

###  3.1 Loss Basin Geometry Comparison

Here, we aim to systematize the observation that the Nuclear estimator can
have slightly higher, but much flatter, minima than the Ridge estimator. To do
so, we simply fix the hyperparameters $\sigma,\lambda$, and consider the
theoretical test error as a function of $\alpha$. We record both the minimal
value over $\alpha$ $m:=min_{\alpha}Err(\alpha)$ as well as the curvature at
the minimum:
$\kappa:=\sqrt{\partial^{2}Err(\alpha)/\partial\alpha^{2}}|_{\alpha=\alpha_{min}}$.
For all cases, we fix $\beta=1$. For the sake of having an interpretable frame
of reference, we report these values as percentage increases relative to the
corresponding values for the Ridge estimator, shown in Figure 3. For the
nuclear case, we can plainly see that the minima of the test loss are nearly
as deep as for the the Ridge (within a few percentage points), however the
minima are often very substantially flatter (corresponding to negative values
in the table). The Spectral estimator tends to also have flatter minima than
Ridge, but with a less favorable depth tradeoff compared to the Nuclear.
Moreover, we see a pronounced effect of $\sigma$, with increasing values
tending to decrease the gap between the depths of Nuclear and Ridge minima,
while also decreasing the difference in curvature. As in Figure 4, we see that
these observations largely carry over to the case of a power-law spectral
density, indicating that they are not simply an incidental property of the
spherical Gaussian model. See A.7 for more technical details on the
calculation of these values.

Figure 3: Quantification of loss basin geometry for the spherical Gaussian
model. Each entry in the table shows the minimal attainable test
error/curvature at the minimum, where each is expressed as a percentage
increase relative to the corresponding value for the Ridge estimator. Figure
4: Same as Figure 3 except for the Diagonal matrix model with power law
spectral density.

###  3.2 Simulation Studies

The purpose of this section is to show that the observations about the loss
basin geometry can have practical consequences in settings where the
regularization strength is selected by cross-validation.

####  3.2.1 Gaussian Predictors

We generated training data by sampling $X_{i}\sim N(0,C(\rho))$, where
$C(\rho):=(1-\rho)I_{d}+\rho\textbf{1}\textbf{1}^{t}$. We then generated a
ground truth coefficient vector by sampling $\beta_{0}\sim N(0,I_{d})$, and
training targets $y_{i}$ by $y_{i}\sim N(\langle
X_{i},\beta_{0}\rangle,\sigma^{2})$. The testing data were generated
similarly, using the same $\beta_{0}$. In all cases, we used $N=100$ training
examples, and $5000$ testing examples. For each combination of hyperparameters
$d,\sigma,\rho$, we generated 100 train/test datasets. For each model
(Spectral, Ridge, Nuclear) and each individual dataset, we used 3-fold cross
validation to select the best-performing $\alpha$ on the training set. We then
refit the model on the entire training set, and evaluated its performance on
the test set. The set of allowable $\alpha$ values considered in the cross-
validation was the same for all models, and consisted of 9 equally
logarithmically spaced values spanning $10^{-4}$ to $10^{6}$.

In Figure 5 we plot the best-performing model for each combination of
hyperparameters, for two different definitions of “best”. In the top row, we
show the model that attains the lowest average test error in each cell. In the
bottom, we show the model that is most likely to “win” on any given dataset.
That is, suppose we evaluate model $m$ on dataset $i$, and obtain a test error
of $MSE_{mi}$. The top row shows $argmin_{m}\mathbb{E}_{j}MSE_{mj}$, while the
bottom row shows $Mode(\\{argmin_{m}MSE_{mj}\\}_{j})$. We see that in terms of
average error, the Ridge and Nuclear are essentially matched, with Spectral
being a distant 3rd. The Nuclear however is overwhelmingly likely to be the
winner on a given dataset.

Figure 5: Best performing models with cross-validation on Gaussian predictors.
Top: best average error. Bottom: highest probability of winning. Each cell is
an aggregate over 100 datasets.

####  3.2.2 Regression of Nonlinear functions using Random Fourier features

In order to generalize our observations beyond the Gaussian features case, we
consider a more complex setup in which the true relationship between the
predictors and the targets is non-linear, and the regression is performed in
the feature space of a universal approximating kernel. To be more specific,
the process for generating data was as follows. We first pick integers
$d,d_{rbf}$ which define the dimension of the observed space, and of the
kernel feature map, respectively. We sampled a training matrix $X$ of size
$N\times d$, in which each entry has variance $1/d_{rbf}$. We then sampled a
non-linear function $f:\mathbb{R}^{d}\to\mathbb{R}$ in a manner described
below, and generated training targets as $y_{i}=f(X_{i})+\sigma\epsilon_{i}$,
where $X_{i}$ is the ith row of $X$ (i.e. ith observation) and $\epsilon_{i}$
is unit Gaussian. Testing data was generated similarly, using the very same
function $f$ (and setting $\sigma=0$ as in the previous analyses). Since the
relation between $X$ and $Y$ is non-linear, we do not directly fit the models
on $X$, but rather first fit a Random Fourier features (RFF) model (Rahimi &
Recht, 2008) on $X$, using $d_{rbf}$ features 666The RFF fit was done using
the class sklearn.kernel_approximation.RBFSampler..We then fit the linear
models on $RFF(X),Y$, where $RFF(X)$ is the $N\times d_{rbf}$ matrix obtained
by applying the Random Fourier features mapping. To estimate test error, we
first transform $X_{test}$ using the same RFF mapping, and consider
$|RFF(X_{test})\hat{\beta}-Y_{test}|^{2}$. Finally, the non-linear function
$f$ was defined as follows:

$f(x)=\sum_{k=1}^{d_{rbf}}cos(2\pi k\langle x,v_{i}\rangle)/k^{2}$ (13)

where $v_{i}$ are sampled uniformly from the unit sphere
$S^{d-1}\subset\mathbb{R}^{d}$. We fixed $d=10$ and $N=100$ in all cases. For
each value of $d_{rbf}$ and $\sigma$, we generated 100 datasets in this way.
We did not consider the Spectral estimator for this experiment, because it
does not naturally accomodate the overparametrized case.As in Figure 6, we see
that the cross-validated Nuclear estimator can again outperform the cross-
validated Ridge, with the effect becoming more pronounced for larger values of
$\sigma$.

Figure 6: Test set performance for cross-validated Nuclear estimator on the
Random Fourier features data. Left shows average test error expressed as a
ratio relative to average Ridge test error. Right shows probability of
“winning” on a given dataset. Each dot represents an aggregate over 100
datasets.

##  4 Related Work

Our work relates to several distinct areas. On the one hand, there are
questions from classical statistics. For instance, several other sorts of
generalized Gauss-Markov theorems have appeared in the literature, although
they generalize along different axes than our formulation. for example
allowing non-spherical errors (Kourouklis & Paige, 1981),predictor
collinearity (Lewis & Odell, 1966), random effects (Shaffer, 1991),or non-
linear estimators (Hansen, 2022) .Similarly,the Nuclear and Spectral
estimators bear a very strong resemblance to various notions from classical
statistics, most notably Principal Components Regression (PCR) (Kendall,
1957), Stein Shrinkage (Stein, 1962) respectively. Finally, the work of
Hocking et al. (1976) is particularly relevant and worthy of further
discussion here-we give a more detailed comparison in A.2.

Lately, there there has also been renewed interest in the study of high-
dimensional regression models and the effects of various kinds of
regularization thereupon. Most such works focus on Ridge or kernel Ridge
regression ((Mei & Montanari, 2020; Dobriban & Wager, 2018; Canatar et al.,
2021; Dicker, 2016; Belkin et al., 2020; Hu & Lu, 2022; Nakkiran et al.,
2021)). Ridge regression has also proven a powerful lens to understand other
techniques such as dropout (Wager et al., 2013), data augmentation(Lin et al.,
2023) and early stopping (Ali et al., 2019). This pervasiveness of the ridge
underscores our claim that a detailed understanding of other kinds of
regularization than L2 in the linear case could be used to study more complex
regularizations in non-linear settings.

Some other work has analyzed other regularizations than ridge such as Bayati &
Montanari (2011); Samet et al. (2013)for Lasso. Theoretical properties of
nuclear norm have mainly been studied in the context of matrix completion, in
which it is often used as a surrogate for matrix rank (Candès & Recht, 2012;
Koltchinskii et al., 2011), but see Hu et al. (2021) for other applications.
Another setting in which rank penalties appear in regression problems is Low
Rank Regression (Bunea et al., 2011; Izenman, 2008), although this is fairly
different than our setting as it has to do with regularizing the outputs of a
regression problem with multiple dependent variables.There is also some modern
work on theoretical properties of the Principal Components Regression
estimator (Xu & Hsu, 2019); optimizing the number of PCR components has a
similar flavor to optimizing $\alpha$ in the Nuclear estimator.

Another line of related work is that while we have here studied the behavior
of flat minima in the context of linear models and cross-validation, this same
concept has also been quite popular as an explanatory tool in the context of
generalization of neural networks (Baldassi et al., 2021; Mulayoff & Michaeli,
2020; Hochreiter & Schmidhuber, 1997), in which the randomness in the
minimization procedure comes from the stochasticity of gradient descent.

##  5 Discussion

We have introduced two simple Linear models:Spectral and Nuclear regression,
and have shown how the forms of both of these models, together with Ridge
regression, may all be derived from an extension of Gauss-Markov Theorem for
Schatten norm constraints on the bias operator. By considering theoretical
error curves as a function of regularization strength, we observed in multiple
different random matrix ensembles that the Nuclear model can often attain
minima that are much flatter than those of the Ridge model, while being nearly
as deep. Finally, we showed that this tradeoff can have practical consequences
in a cross-validation setting, for both Gaussian and non-Gaussian features.

It should be noted that the particulars of the cross-validation results depend
on the number of $\alpha$ values used in the grid search, among other factors;
and as noted before, the effect of the curvature at the minimum decreases as
more values of $\alpha$ are used. However, increasing the number of $\alpha$
values has other potential costs such as risk of overfitting or increased
computational cost. But more fundamentally,our intention with the simulation
results is not to claim that the Nuclear estimator should always be preferred
over the Ridge, but rather to show that it can outperform in some practically
relevant setups, and to use our theoretical analysis to give a precise
characterization of why this happens.

There are many potential future directions of this work, such as further
understanding of the Spectral and Nuclear regularizations in the kernelized or
overparametrized regime, or using insights from the Nuclear case to understand
other situations that involve rank constraints such as matrix completion. More
broadly, we hope that this framework will provide a test case to understand
less common regularizations such as $\|\cdot\|_{1}$ and $\|\cdot\|_{\infty}$,
with the goal of transferring the insights gained thereby to the study of
high-dimensional non-linear models.

## References

  * Ali et al. (2019) A Ali, JZ Kolter, and RJ Tibshirani.  A continuous-time view of early stopping for least squares.  _AISTATS_ , 2019.
  * Appell (1925) P Appell.  Sur les fonctions hypergéométriques de plusieurs variables.  _Mémoir. Sci. Math. Paris: Gauthier-Villars_ , 1925.
  * Bai & Silverstein (2010) Z Bai and J Silverstein.  _Spectral Analysis of Large Dimensional Random Matrices_.  Springer, 2010.
  * Bailey (1934) WN Bailey.  On the reducibility of appell’s function $f_{4}$.  _Quart. J. Math_ , 1934.
  * Baldassi et al. (2021) C Baldassi, C Lauditi, EM Malatesta, G Perugini, and R Zecchina.  Unveiling the structure of wide flat minima in neural networks.  _Phys Rev Letters_ , 2021.
  * Bayati & Montanari (2011) M Bayati and A Montanari.  The dynamics of message passing on dense graphs, with applications to compressed sensing.  _IEEE Transcations on Information Theory_ , 2011.
  * Belkin et al. (2019) M Belkin, D Hsu, S Ma, and S Mandal.  Reconciling modern machine learning practice and the bias-variance trade-off.  _Proc Nat Aca Sci_ , 2019.
  * Belkin et al. (2020) M Belkin, D Hsu, and J Xu.  Two models of double descent for weak features.  _preprint arXiv:1903.07571_ , 2020.
  * Bunea et al. (2011) F Bunea, Y She, and M Wegkamp.  Optimal selection of reduced rank estimators of high-dimensional matrices.  _Ann. Statist._ , 2011.
  * Canatar et al. (2021) A Canatar, B Bordelon, and C Pehlevan.  Spectral bias and task-model alignment explain generalization in kernel regression and infinitely wide neural networks.  _Nature communications_ , 2021.
  * Candès & Recht (2012) E Candès and B Recht.  Exact matrix completion via convex optimization.  _Communications of the ACM_ , 2012.
  * Dicker (2016) LH Dicker.  Ridge regression and asymptotic minimax estimation over spheres of growing dimension.  _International Statistical Institute and the Bernoulli Society for Mathematical Statistics and Probability_ , 2016.
  * Dobriban & Wager (2018) E Dobriban and S Wager.  High-dimensional asymptotics of prediction: Ridge regression and classification.  _Annals of statistics_ , 2018.
  * Fan (1951) Ky Fan.  Maximum properties and inequalities for the eigenvalues of completely continuous operators.  _Proc. Nat. Acad. Sci. U.S.A._ , 1951.
  * Goodfellow et al. (2016) I Goodfellow, Y Bengio, and A Courville.  _Deep learning_.  MIT press, 2016.
  * Hansen (2022) BE Hansen.  A modern gauss–markov theorem.  _Econometrica_ , 2022.
  * Hastie et al. (2019) T Hastie, A Montanari, S Rosset, and RJ Tibshirani.  Surprises in high-dimensional ridgeless least squares interpolation.  _arxiv:abs/1903.08560_ , 2019.
  * Hochreiter & Schmidhuber (1997) S Hochreiter and J Schmidhuber.  Flat minima.  _Neural Computation_ , 1997.
  * Hocking et al. (1976) R. R. Hocking, F. M. Speed, and M. J. Lynn.  A class of biased estimators in linear regression.  _Technometrics_ , 1976.
  * Hu & Lu (2022) H Hu and YM Lu.  Sharp asymptotics of kernel ridge regression beyond the linear regime.  _arXiv:2205.06798_ , 2022.
  * Hu et al. (2021) Z Hu, F Nie, R Wang, and X Li.  Low rank regularization: A review.  _Neural Networks_ , 2021.
  * Izenman (2008) AJ Izenman.  Modern multivariate statistical techniques: Regres- sion, classification, and manifold learning.  _springer_ , 2008.
  * Jacot et al. (2020) A Jacot, B Simsek, F Spadaro, C Hongler, and F Gabriel.  Kernel alignment risk estimator: risk prediction from training data.  _In Advances in Neural Information Processing Systems_ , 2020.
  * Johansson et al. (2022) A Johansson, N Engsner, C Strannegard, and P Mostad.  Exact spectral norm regularization for neural networks.  _arxiv:2206.13581.pdf_ , 2022.
  * Kamalakara et al. (2022) S Rao Kamalakara, A Locatelli, B Venkitesh, J Ba, Y Gal, and AN Gomez.  Exploring low rank training of deep neural networks.  _arXiv:2209.13569_ , 2022.
  * Kendall (1957) MG Kendall.  _A Course in Multivariate Analysis_.  London: Charles Griffin, 1957.
  * Koltchinskii et al. (2011) V Koltchinskii, K Lounici, and AB Tsybakov.  Nuclear-norm penalization and optimal rates for noisy low-rank matrix completion.  _Ann. Statist_ , 2011.
  * Kourouklis & Paige (1981) S Kourouklis and CC Paige.  A constrained least squares approach to the general gauss- markov linear model.  _Journal of American Statistical Assoc._ , 1981.
  * Lewis & Odell (1966) TO Lewis and PL Odell.  A generalization of the gauss-markov theorem.  _Journal of American statistical assoc._ , 1966.
  * Lin et al. (2023) CH Lin, C Kaushik, EL Dyer, and V Muthukumar.  The good, the bad and the ugly sides of data augmentation: An implicit spectral regularization perspective.  _arXiv:2210.05021_ , 2023.
  * Liu et al. (1997) Y Liu, P Cizeau, M Meyer, CK Peng, and HE Stanley.  Correlations in economic time series.  _Physica A: Statistical Mechanics and its Applications_ , 1997.
  * Mei & Montanari (2020) S Mei and A Montanari.  The generalization error of random features regression: Precise asymptotics and double descent curve.  _arxiv.org/pdf/1908.05355.pdf_ , 2020.
  * Mianjy & Arora (2019) P Mianjy and R Arora.  On dropout and nuclear norm regularization.  _International conference on machine Learning_ , 2019.
  * Mulayoff & Michaeli (2020) R Mulayoff and T Michaeli.  Unique properties of flat minima in deep networks.  _Proceedings of the 37th International Conference on Machine Learning_ , 2020.
  * Nakkiran et al. (2021) P Nakkiran, P Venkat, SM Kakade, and T Ma.  Optimal regularization can mitigate double descent.  _In International Conference on Learning Representations_ , 2021.
  * Qin & Colwell (2018) C Qin and LJ Colwell.  Power law tails in phylogenetic systems.  _Proc. Nat. Aca. Sci_ , 2018.
  * Rahimi & Recht (2008) A Rahimi and B Recht.  Weighted sums of random kitchen sinks: Replacing minimization with randomization in learning.  In _NeurIPS_ , 2008.
  * Ruderman & Bialek (1994) DL Ruderman and W Bialek.  Statistics of natural images: Scaling in the woods.  _Physical review letters_ , 1994.
  * Samet et al. (2013) O Samet, C Thrampoulidis, and B Hassibi.  The squared-error of generalized lasso: A precise analysis.  _Allerton Conference on Communication, Control, and Computing_ , 2013.
  * Saxe et al. (2013) AM Saxe, JL McClelland, and S Ganguli.  Exact solutions to the nonlinear dynamics of learning in deep linear neural networks.  _arXiv:1312.6120_ , 2013.
  * Shaffer (1991) JP Shaffer.  The gauss—markov theorem and random regressors.  _The American Statistician_ , 1991.
  * Srivastava et al. (2014) N Srivastava, G Hinton, A Krizhevsky, I Sutskever, and R Salakhutdinov.  Dropout: A simple way to prevent neural networks from overfitting.  _The Journal of Machine Learning Research_ , 2014.
  * Stein (1962) J Stein.  Multiple regression.  _Contributions to Probability and Statistics. Essays in Honor of Harold Hotelling_ , pp. 424–443, 1962.
  * Stringer et al. (2019) C Stringer, M Pachitariu, N Steinmetz, M Carandini, and KD Harris.  High-dimensional geometry of population responses in visual cortex.  _Nature_ , 2019.
  * Wager et al. (2013) S Wager, S Wang, and P Liang.  Dropout training as adaptive regularization.  _arXiv:1307.1493v2_ , 2013.
  * Wooldride (2015) JM Wooldride.  _Introductory Econometrics: A Modern Approach_.  Cengage Learning, 2015.
  * Xu & Hsu (2019) J Xu and D Hsu.  On the number of variables to use in principal component regression.  _Neurips_ , 2019.
  * Yang et al. (2020) H Yang, M Tang, W Wen, F Yan, D Hu, A Li, H Li, and Y Chen.  Learning low-rank deep neural networks via singular vector orthogonality regularization and singular value sparsification.  _arXiv:2004.09031_ , 2020.
  * Yoshida & Miyato (2017) Y Yoshida and T Miyato.  Spectral norm regularization for improving the generalizability of deep learning.  _arXiv:1705.10941_ , 2017.

##  Appendix A appendix

###  A.1 Code availability

Code for the simulations is available at
https://github.com/SimonSegert/specreg.

###  A.2 Further comparison with Hocking et. al.

The paper, similarly to us in spirit, introduces a class of biased linear
estimators defined by a family of constrained optimization problems, and
derive in this way the forms of Ridge, Stein Shrinkage, and PCR regressors.
However, a closer inspection reveals that their setup is considerably
different from ours.

For simplicity, we will make the assumption that $G=X^{T}X$ is diagonal, as
this is sufficient to illustrate the differences. The form of estimators
considered in Hocking et. al. is

$\hat{\beta}=v\circ\Sigma$ (14)

where as before, $\Sigma$ is the vector of eigenvalues of $G$ and $\circ$ is
Hadamard (element-wise) product. Immediately we see a difference that their
setup has only $d$ free variables, from the vector $p$, whereas we have
$N\times p$ free variables in the matrix $L$.

By a bit of algebra,

$\displaystyle\hat{\beta}$ $\displaystyle=$ $\displaystyle\Sigma\circ
v=(\Sigma\circ v\oslash X^{T}Y)\circ X^{T}Y$ (15) $\displaystyle=$
$\displaystyle diag(\Sigma)diag(v)diag(X^{T}Y)^{-1}X^{T}Y$ (16)
$\displaystyle:=$ $\displaystyle DX^{T}Y$ (17)

where $\oslash$ is element-wise division, and $D$ is a diagonal matrix. Thus
they effectively assume from the outset that $\hat{\beta}=DX^{T}Y$ for some
diagonal matrix $D$, whereas we a priori allow $\hat{\beta}=LY$ for
$L\in\mathbb{R}^{N\times d}$ and derive that $L$ is necessarily of this form.

They then show that each of the aforementioned estimators can be derived as
the solution to an optimization of the form $\min_{v\in C}L(v)$, where $L$ is
a loss function and $C$ is a constraint set, that depends on the estimator.
The forms of the constraint sets are rather ad-hoc and the authors do not
appear justify the form from more basic principles, or to relate them to each
other in a meaningful way. By contrast, the forms of our constraints arise
naturally from a standard family of matrix norms.

Finally, they do not precisely characterize the bias of each of the
estimators, whereas our results characterize the estimators as being optimal
for a fixed level of bias.

###  A.3 Rule of Thumb for depth vs. curvature

In order to get insight into the tradeoff between curvature and depth, we
consider a highly simplified setup which nonetheless captures some key
features of selecting $\alpha$ through cross-validation. Let us model the
function of average test error vs regularization strength as a parabola
$f(x)=\kappa^{2}x^{2}/2+\mu$, for $-\delta\leq x\leq\delta$. Thus $\kappa$
controls the curvature and $\mu$ controls the depth of the minimum. We will
model cross-validation as taking n iid samples $X_{i}\sim
Unif(-\delta,\delta)$ and then estimating the minimum of $f$ as
$\hat{\mu}=min_{i}(f(X_{i}))$. Note that in a real cross-validation setup
there is the additional complication that we cannot perfectly observe the
value $f(x)$ due to finite-sample variability; we do not model this effect
here.

What is the distribution of $\hat{\mu}$? Well, for
$\mu<y<\delta^{2}\kappa^{2}/2+\mu=f_{max}$, we have

$\mathbb{P}(f(X_{i})>y)=\mathbb{P}(|X_{i}|>\sqrt{2(y-\mu)}/\kappa={\frac{1}{2\delta}}(2\delta-2\sqrt{2(y-\mu)}/\kappa)$
(18)

So,

$\mathbb{P}(\hat{\mu}>y)=\mathbb{P}(f(X_{i})>y)^{n}=(1-\delta^{-1}\sqrt{2(y-\mu)}/\kappa)^{n}$
(19)

The expected value is

$\mathbb{E}\hat{\mu}=\int_{0}^{\infty}\mathbb{P}(\mu>y)dy=\mu+\int_{\mu}^{\delta^{2}\kappa^{2}/2+\mu}(1-\sqrt{(2(y-\mu)\kappa^{-2}\delta^{-2}})^{n}dy$
(20)

Letting $u=\sqrt{2(y-\mu)/(\kappa^{2}\delta^{2})}$,the integral becomes

$\mu+\kappa^{2}\delta^{2}\int_{0}^{1}(1-u)^{n}udu=\mu+\kappa^{2}\delta^{2}Beta(2,n+1)=\mu+{\frac{(\kappa\delta)^{2}}{(n+1)(n+2)}}$
(21)

In the large sample limit, the contribution of the curvature vanishes, and we
can exactly find the minimum. However, for finite samples this is not the
case, and it could be that for two different parabolas the true minima satisfy
$\mu_{1}<\mu_{2}$ while the estimated minima satisfy
$\mathbb{E}\hat{\mu_{1}}>\mathbb{E}\hat{\mu_{2}}$.

###  A.4 proof of main theorem

Let us first make the simple but important observation that a minimizer in
Equation 3 actually exists (since the constraint set is non-compact this is
not completely immediate). One way to see this is to note that the problem is
of the form $\min_{x\in D}|x|^{2}$ for some closed set $D$, and problems of
this form always have a minimizer (even if $D$ is non-compact). Combined with
the observation that the objective function is strictly convex, we conclude
that the problem in Equation 3 has exactly one  minimizer.

The basic strategy of the proof is now to derive certain properties of the
minimizer, and add these properties as further constraints until we get
something tractable.

######  Lemma 3.

The minimizer $L_{opt}$ of the bias-constrained problem (Equation 3) takes the
form $QX^{T}$ for some $d\times d$ matrix $Q$

###### Proof.

Let $X_{\perp}$ be any $N\times(N-d)$ matrix whose columns for a basis for the
orthogonal complement of $Colspace(X)$. We may assume that the columns are
orthonormal; $X_{\perp}^{T}X_{\perp}=I_{N-d}$. Note that the concatenated
matrix $[X,X_{\perp}]$ is invertible. Therefore, we can express the optimum as
$L=M\left(\begin{array}[]{c}X^{T}\\\ X_{\perp}^{T}\end{array}\right)$ for some
$d\times N$ matrix $M$. Writing $M$ in block form $M=[Q,Q_{\perp}]$ where
$Q\in\mathbb{R}^{d\times d}$ and $Q_{\perp}\in\mathbb{R}^{d\times(N-d)}$, we
have $L=QX^{T}+Q_{\perp}X_{\perp}^{T}$. Note that $X^{T}X_{\perp}=0$;
therefore

$\displaystyle\|LL^{T}\|_{F}^{2}$ $\displaystyle=$
$\displaystyle\|QX^{T}\|_{F}^{2}+\|Q_{\perp}X_{\perp}^{T}\|_{F}^{2}$ (22)
$\displaystyle=$
$\displaystyle\|QX^{T}\|_{F}^{2}+Tr(Q_{\perp}X_{\perp}^{T}X_{\perp}Q_{\perp}^{T})$
(23) $\displaystyle=$
$\displaystyle\|QX^{T}\|_{F}^{2}+Tr(Q_{\perp}Q_{\perp}^{T})$ (24)
$\displaystyle=$ $\displaystyle\|QX^{T}\|_{F}^{2}+\|Q_{\perp}\|_{F}^{2}$ (25)

where $\|\cdot\|_{F}$ is the Frobenius norm. Similarly, the bias
$LX-I=QX^{T}X-I$ does not depend on $Q_{\perp}$. Thus, any non-zero value of
$Q_{\perp}$ will strictly increase the value of the objective relative to
setting $Q_{\perp}=0$, without having any effect on the constraint. ∎

######  Corollary 3.1.

The optimal solution to Equation 3 takes the form $L=(I+Q)G^{-1}X^{T}$ where

$Q=argmin_{Q^{\prime}\in\mathbb{R}^{d\times d}}Tr(Q^{\prime}G^{-1}Q^{\prime
T})/2+Tr(Q^{\prime}G^{-1})+\alpha^{-1}\|Q^{\prime}\|_{p}$ (26)

and $\alpha\geq 0$ is determined by $C$.

###### Proof.

By the Lemma, it is no loss of generality to assume that $L$ has the indicated
form. Now plug in to Equation 3 and simplify. The $\alpha$ term is just
converting the constraint to a Lagrange multiplier. ∎

By the discussion above, we conclude that there is exactly one minimum of
Equation 26.

######  Proposition 3.1.

If $Q$ is the solution to the equation 26, then $Q$ is symmetric and commutes
with $G$

Before giving the proof, we first present a more technical matrix lemma.

######  Lemma 4.

For any square matrix $M$ and $p\geq 1$, $\|M_{d}\|_{p}\leq\|M\|_{p}$, where
$M_{d}$ is the diagonal part (i.e., the matrix with all non-diagonal entries
set to zero).

###### Proof.

Evidently the singular values of $M_{d}$ coincide with the absolute values of
the diagonal entries. By an inequality of Ky Fan (Fan, 1951),

$\sum_{i\leq k}\sigma_{i}(M_{d})\leq\sum_{i\leq k}\sigma_{i}(M)$ (27)

for any $1\leq k\leq d$, where $\sigma_{i}$ denotes the singular values,
ordered from largest to smallest. The lemma now follows from Schur convexity
of the Euclidean p-norm. ∎

###### Proof.

(of proposition 3.1) Let $G^{-1}=UDU^{T}$ be the diagonalization, where $U$ is
an orthogonal matrix and $D$ is diagonal. The objective is

$\displaystyle Q$ $\displaystyle=$ $\displaystyle
argmin_{Q}Tr(Q^{\prime}UDU^{T}Q^{\prime
T})/2+Tr(Q^{\prime}UDU^{T})+\alpha^{-1}\|Q^{\prime}\|_{p}$ (28)
$\displaystyle=$ $\displaystyle argmin_{Q}Tr(U^{T}Q^{\prime}UDU^{T}Q^{\prime
T}U)/2+Tr(U^{T}Q^{\prime}UD)+\alpha^{-1}\|U^{T}Q^{\prime}U\|_{p}$ (29)
$\displaystyle UQU^{T}$ $\displaystyle=$ $\displaystyle
argmin_{P}Tr(PDP^{T})/2+Tr(PD)+\alpha^{-1}\|P\|_{p}$ (30)

where we reparametrized as $P:=U^{T}Q^{\prime}U$. Evidently, the proposition
will follow if we can show that the minimal $P$ is diagonal. Let
$P=P_{d}+P_{od}$ be the decomposition into diagonal and off-diagonal
parts777By definition an off-diagonal matrix is one with all zeros along the
diagonal.. Plugging into the objective

$\displaystyle Obj(P)$ $\displaystyle=$ $\displaystyle
Tr((P_{d}+P_{od})D(P_{d}+P_{od}^{T}))/2+Tr(P_{d}D)+Tr(P_{od}D)+\alpha^{-1}\|P\|_{p}$
(32) $\displaystyle=$ $\displaystyle
Tr(P_{d}DP_{d})/2+Tr(P_{d}D)+Tr(P_{od}DP_{od}^{T})/2+Tr(P_{od}D)+\alpha^{-1}\|P\|_{p}$
(33) $\displaystyle+$ $\displaystyle Tr(P_{d}DP_{od}^{T})$ (34)

Now, we note that in general if $A$ is diagonal, and $B$ is off-diagonal,then
$Tr(AB)=0$. So the expression simplifies to

$Tr(P_{d}DP_{d})/2+Tr(P_{d}D)+Tr(P_{od}DP_{od}^{T})/2+\alpha^{-1}\|P\|_{p}$
(35)

Now, the third term is non-negative because it is the trace of a PSD matrix;
thus $Obj(P)\geq Tr(P_{d}DP_{d})/2+Tr(P_{d}D)+\alpha^{-1}\|P\|_{p}$. By Lemma
4, $\|P\|_{p}\geq\|P_{d}\|_{p}$, and therefore

$Obj(P)\geq Obj(P_{d})$ (36)

However, $P$ was assumed to be the minimum, which implies by strict convexity
that $P=P_{d}$ ∎

By the above, we now know that the optimal $Q$ must satisfy $Q=Udiag(q)U^{T}$,
where $q$ denotes the vector of eigenvalues and $U$ is the matrix of
eigenvectors of $G^{-1}$. Since $\hat{G}^{-1}=(1+Q)G^{-1}$ in the notation of
the theorem statement, we have shown the first two claims, namely that
$\hat{G}$ is symmetric and commutes with $G$.

By plugging into the objective in 26 and simplifying, we see that

$q=argmin_{x\in\mathbb{R}^{d}}{\frac{1}{2}}\sum_{i}x_{i}^{2}/\sigma_{i}^{2}+\sum_{i}x_{i}/\sigma_{i}^{2}+\alpha^{-1}|x|_{p}$
(37)

where $\sigma_{i}^{2}$ are the eigenvalues of $G$, and $|\cdot|_{p}$ denotes
the Euclidean p norm of a vector.

To show the last claim that $\hat{G}-G$ is non-negative definite, it is enough
to show that $0\geq q_{i}\geq-1$,since
$\lambda_{i}(\hat{G})=\lambda_{i}(G)/(1+q_{i})$.It is easy to see that the
eigenvectors of $Q$ must be non-positive (since if any eigenvector is
positive, then switching the sign will leave the first and third terms alone,
while strictly decreasing the second term). To see the second inequality, we
first add a suitable constant to the objective and rewrite it as
${\frac{1}{2}}\sum_{i}(x_{i}+1)^{2}/\sigma_{i}^{2}+\alpha^{-1}|x|_{p}$.
Supposing that $q_{i}<-1$ for some $i$, let us replace it with
$q_{i}^{\prime}=-2-q_{i}$. Doing so does not change the first term (i.e.
$(q_{i}+1)^{2}=(q_{i}^{\prime}+1)^{2}$), however it strictly decreases the
second term, since $|q_{i}|^{\prime}<|q_{i}|$ 888e.g., square both sides,
contradicting the minimality of $q$999This argument doesn’t quite go through
if $p=\infty$, since in that case decreasing the magnitude of a component of
$q$ might not decrease the norm. But that’s fine because we will derive the
exact solution for $p=\infty$ shortly.

Until now, we have not made any assumption about $p$ except that $p\geq 1$. At
this point, we separately analyze each of the three special cases
$p=1,2,\infty$.

Before doing so, however, we first present the following well-known and
elementary calculation, which we will employ several times in what follows:

######  Proposition 4.1.

Let $y,\tau>0$ then

$argmin_{x\in\mathbb{R}}{\frac{1}{2}}(x+y)^{2}+\tau|x|=min(\tau-y,0)$ (38)

###### Proof.

It is clear that the optimal $x$ cannot be positive, so we want to compute

$argmin_{x\leq 0}{\frac{1}{2}}(x+y)^{2}-\tau x=argmin_{x\leq
0}{\frac{1}{2}}(x+y-\tau)^{2}+{\frac{y^{2}-(y-\tau)^{2}}{2}}$ (39)

where we completed the square. The second term does not depend on $x$ and
therefore has no effect on the argmin. Now the formula is immediate. If
$\tau-y<0$, then the minimum is plainly attained at $x=\tau-y<0$. If
$\tau-y>0$, then the parabola is monotonically decreasing on the interval
$(-\infty,0)$, so the minimum is attained at $x=0$.

∎

Nuclear case (p=1)

The objective 37 splits into a sum of separable one-dimensional problems:

$argmin_{x\in\mathbb{R}}{\frac{1}{2}}x^{2}/\sigma_{i}^{2}+x/\sigma_{i}^{2}+\alpha^{-1}|x|$
(40)
$argmin_{x\in\mathbb{R}}{\frac{1}{2}}(x+1)^{2}+\sigma_{i}^{2}\alpha^{-1}|x|$
(41)

By Proposition4.1,

$q_{i}=min(\sigma_{i}^{2}/\alpha-1,0)$ (42)

The eigenvalues of $\hat{G}$ are thus given by

$\lambda_{i}(\hat{G})={\frac{\sigma_{i}^{2}}{1+q_{i}}}={\frac{\sigma_{i}^{2}}{min(\sigma_{i}^{2}/\alpha,1)}}=max(\sigma_{i}^{2},\alpha)$
(43)

as claimed.

Frobenius case (p=2)

This is a simple calculus exercise and omitted.

Spectral case (p=$\infty$)

W know that there exists exactly one minimizer $q_{opt}\in\mathbb{R}^{d}$ of
37. Fix some $C>max(|q_{opt}|_{\infty},max_{i}(1+\sigma_{i}^{2}/\alpha))$ and
consider the problem

$\min_{q:|q|_{\infty}\leq
C}{\frac{1}{2}}\sum_{i}q_{i}^{2}/\sigma_{i}^{2}+\sum_{i}q_{i}/\sigma_{i}^{2}+\alpha^{-1}|q|_{\infty}$
(44)

This clearly has the same minimum as the original unconstrained problem.

Note the following identity:

$\alpha^{-1}|q|_{\infty}=max_{y:|y|_{1}\leq\alpha^{-1}}\langle q,y\rangle$
(45)

By Von Neumann’s minimax theorem, we can interchange the min and the max. So
the optimal value of the objective is:

$max_{y:|y|_{1}\leq\alpha^{-1}}min_{q:|q|_{\infty}\leq
C}{\frac{1}{2}}\sum_{i}q_{i}^{2}/\sigma_{i}^{2}+\sum_{i}q_{i}/\sigma_{i}^{2}+\langle
q,y\rangle$ (46)

Now, the set $\\{q:|q|_{\infty}\leq C\\}$ is geometrically a product of
intervals $[-C,C]^{d}$. Thus, we see that the inner objective splits into a
sum of uncoupled 1-dimensional problems:

$q_{i}=argmin_{q_{i}\in[-C,C]}{\frac{1}{2}}q_{i}^{2}/\sigma_{i}^{2}+q_{i}/\sigma_{i}^{2}+q_{i}y_{i}=argmin_{q_{i}\in[-C,C]}{\frac{1}{2\sigma_{i}^{2}}}(q_{i}+1+\sigma_{i}^{2}y_{i})^{2}-{\frac{1}{2\sigma_{i}^{2}}}(1+\sigma_{i}^{2}y_{i})^{2}$
(47)

By assumption on $C$, $1+\sigma_{i}^{2}y_{i}\leq C$, and therefore for any
fixed $y_{i}$, the minimum of the inner problem is attained at

$q_{i}=-1-\sigma_{i}^{2}y_{i}$ (48)

, with minimal value equal to
$-{\frac{1}{2\sigma_{i}^{2}}}(1+\sigma_{i}^{2}y_{i})^{2}=-{\frac{\sigma_{i}^{2}}{2}}(\sigma_{i}^{-2}+y_{i})^{2}$.
Plugging back in to 46, we need to solve

$min_{y:|y|_{1}\leq\alpha^{-1}}{\frac{1}{2}}\sum_{i}\sigma_{i}^{2}(\sigma_{i}^{-2}+y_{i})^{2}$
(49)

For this, we introduce a Lagrange multiplier $\lambda$, upon which the
objective splits again into a sum of uncoupled 1d problems:
$y_{i}=argmin_{y\in\mathbb{R}}{\frac{1}{2}}\sigma_{i}^{2}(\sigma_{i}^{-2}+y)^{2}+\lambda|y|$.

Using Proposition 4.1, we derive the solution

$y_{i}=\sigma_{i}^{-2}min(\lambda-1,0)=\sigma_{i}^{-2}(min(\lambda,1)-1)$ (50)

where $\lambda$ is chosen to satisfy the original constraint
$\sum_{i}|y_{i}|\leq\alpha^{-1}$. Using the relation 48 between $q_{i}$ and
$y_{i}$

$q_{i}=-1-(min(\lambda,1)-1)=-min(\lambda,1)$ (51)

for the solution to the problem 44. Since we took $C$ to be large enough to
contain the solution to the unconstrained problem, we conclude that this is
also the solution to the unconstrained problem (i.e. $C=\infty$). Since
$x_{i}$ are the eigenvalues of $Q$, we conclude that the eigenvalues of
$\hat{G}$ are

$\lambda_{i}(\hat{G})={\frac{\sigma_{i}^{2}}{1+q_{i}}}={\frac{\sigma_{i}^{2}}{1-min(\lambda,1)}}$
(52)

Clearly the denominator lies in $[0,1]$, therefore we recover the claimed form
in which all eigenvalues of $\hat{G}$ are obtained by scalar multiplication
with some factor $>1$. Note that the case $\lambda>1$ corresponds to
multiplication by infinity, i.e. setting $\hat{G}^{-1}$ (and thus
$\hat{\beta})$ to zero.

###  A.5 Proof of 2.3

We first give the definition of the Appel hypergeometric function $F_{1}$ for
reference.

The function is typically defined as

$F_{1}(\alpha,\beta,\beta^{\prime},\gamma,x,y)=\sum_{m=0}^{\infty}\sum_{n=0}^{\infty}{\frac{(\alpha)_{m+n}(\beta)_{m}(\beta^{\prime})_{n}}{m!n!(\gamma)_{m+n}}}x^{m}y^{n}$
(53)

Here $(\cdot)_{m}$ is the Pochammer symbol. The series is absolutely
convergent for $|x|,|y|<1$, and arbitrary
$\alpha,\beta,\beta^{\prime},\gamma$. In 2.3 we may possibly need to evaluate
it at some $x$ outside of this range; we do this by appropriate analytic
continuation, discussed further below.

To prove the formula, we note that the extremal cases follow straightforwardly
from the integral formula. So in what follows we will assume that
$\alpha\in[\lambda_{-},\lambda_{+}]$.

We rewrite the integral

$\displaystyle{\frac{Err}{\lambda}}-\beta^{2}$ $\displaystyle=$
$\displaystyle\int_{\lambda_{-}}^{\lambda_{+}}\beta^{2}x^{2}f_{\alpha}(x)^{-2}-2\beta^{2}xf_{\alpha}^{-1}(x)+\sigma^{2}xf^{-2}_{\alpha}(x)d\mu$
(54) $\displaystyle=$
$\displaystyle\int_{\lambda_{-}}^{\alpha}\beta^{2}\alpha^{-2}x^{2}-2\beta^{2}\alpha^{-1}x+\sigma^{2}\alpha^{-2}xd\mu$
(55) $\displaystyle+$
$\displaystyle\int_{\alpha}^{\lambda_{+}}\beta^{2}-2\beta^{2}+\sigma^{2}x^{-1}d\mu$
(56) $\displaystyle=$
$\displaystyle\beta^{2}\alpha^{-2}\int_{\lambda_{-}}^{\alpha}x^{2}d\mu+(\sigma^{2}\alpha^{-2}-2\beta^{2}\alpha^{-1})\int_{\lambda_{-}}^{\alpha}xd\mu$
(57) $\displaystyle-$
$\displaystyle\beta^{2}(1-F_{\lambda}(\alpha))+\sigma^{2}({\frac{1}{1-\lambda}}-\int_{\lambda_{-}}^{\alpha}x^{-1}d\mu)$
(58) $\displaystyle=$ $\displaystyle
c+{\frac{\beta^{2}}{\alpha^{2}}}I(2,\alpha)+(\sigma^{2}\alpha^{-2}-2\beta^{2}\alpha^{-1})I(1,\alpha)+\beta^{2}F_{\lambda}(\alpha)-\sigma^{2}I(-1,\alpha)$
(59)

where we defined $I(r,\alpha)=\int_{\lambda_{-}}^{\alpha}x^{r}d\mu$.

So we have reduced the theorem to just evaluating $I(r,\alpha)$ for
$r\in\\{-1,1,2\\}$. Now, we plug in the form of the MP density, and make the
variable substitution $u={\frac{x-\lambda_{-}}{\alpha-\lambda_{-}}}$. Clearly
then $u(\alpha-\lambda_{-})+\lambda_{-}=x$.

$\displaystyle 2\pi I(r,\alpha)$ $\displaystyle=$
$\displaystyle\int_{\lambda_{-}}^{\alpha}x^{r-1}\sqrt{(\lambda_{+}-x)(x-\lambda_{-})}dx$
$\displaystyle=$
$\displaystyle\int_{0}^{1}(u(\alpha-\lambda_{-})+\lambda_{-})^{r-1}\sqrt{\lambda_{+}-\lambda_{-}-u(\alpha-\lambda_{-})}\sqrt{u(\alpha-\lambda_{-})}(\alpha-\lambda_{-})du$
$\displaystyle=$
$\displaystyle(\alpha-\lambda_{-})^{3/2}\lambda_{-}^{r-1}\sqrt{\lambda_{+}-\lambda_{-}}\int_{0}^{1}\sqrt{u}(1+u{\frac{\alpha-\lambda_{-}}{\lambda_{-}}})^{r-1}\sqrt{1-u{\frac{\alpha-\lambda_{-}}{\lambda_{+}-\lambda_{-}}}}du$

We can now express this in the form given in the proposition by means of the
following formula (Bailey, 1934):

$F_{1}(\alpha,\beta,\beta^{\prime},\gamma,x,y)={\frac{\Gamma(\gamma)}{\Gamma(\alpha)\Gamma(\gamma-\alpha)}}\int_{0}^{1}u^{\alpha-1}(1-u)^{\gamma-\alpha-1}(1-ux)^{-\beta}(1-uy)^{-\beta^{\prime}}du$
(60)

if $\alpha,\gamma-\alpha>0$. We use this formula to define the analytic
continuation in case $|x|>1$ as can happen in formula 2.3. ∎

###  A.6 Proof of Generalization Formulas 2.1 and 2.4

We first consider 2.4. Given the setup in the main text, consider computing
the test error for a fixed training/testing pair:

$|X_{test}(\hat{G}^{-1}X_{tr}^{T})(X_{tr}\beta_{0}+\epsilon)-X_{test}\beta_{0}|^{2}=|X_{test}(\hat{G}^{-1}G-I)\beta_{0}+X_{test}\hat{G}^{-1}X_{tr}^{T}\epsilon|^{2}$
(61) $=|X_{2}DSB\beta_{0}+X_{2}DS\hat{G}^{-1}X_{tr}^{T}\epsilon|^{2}$ (62)

where $S$ is the diagonal matrix containing the noise values $\sqrt{s_{i}}$,
and $D$ is the diagonal matrix containing $\sqrt{\lambda_{i}}$. By rotational
symmetry of the Frobenius norm, when we marginalize $X_{2}$ we get

$|DSB\beta_{0}+DS\hat{G}^{-1}X_{tr}^{T}\epsilon|^{2}$ (63)

And marginalizing over $\epsilon$ further yields

$|DSB\beta_{0}|^{2}+\sigma^{2}Tr(DS\hat{G}^{-1}G\hat{G}^{-1}SD)$ (64)

where the cross term vanishes because $\mathbb{E}\epsilon=0$. Remembering all
matrices are actually diagonal, we obtain

$n*err=\sum_{i=1}^{d}s_{i}\lambda_{i}(\beta_{0})_{i}^{2}(\lambda_{i}/f_{\alpha}(\lambda_{i})-1)^{2}+\sigma^{2}\sum_{i}s_{i}\lambda_{i}(\lambda_{i}/f_{\alpha}(\lambda_{i})^{2})$
(65)

Since $s_{i}$ are independent of $\lambda_{i}$ and have expectation 1, we can
easily marginalize out, obtaining:

$\sum_{i=1}^{d}\lambda_{i}(\beta_{0})_{i}^{2}(\lambda_{i}/f_{\alpha}(\lambda_{i})-1)^{2}+\sigma^{2}\sum_{i}\lambda_{i}(\lambda_{i}/f_{\alpha}(\lambda_{i})^{2})$
(66)

At this point, the only variable that remains to marginalize over is
$\lambda$. The result is:

$d^{-1}|\beta_{0}|^{2}\sum_{i}E_{\lambda_{1},\dots,\lambda_{n}}\lambda_{i}\left((\lambda_{i}/f_{\alpha}(\lambda_{i})-1)\right)^{2}+\sigma^{2}E_{\lambda_{1},\dots,\lambda_{n}}\sum_{i}\lambda_{i}^{2}/f_{\alpha}(\lambda_{i})^{2}$
(67)

where we used exchangeability of $\lambda_{i}$ to bring out the factor of
$|\beta_{0}|^{2}$. Now the proposition follows by taking the limit
$n\to\infty$ and applying the law of large numbers. ∎

The proof of 2.1 is very similar - as above, the idea is to use rotational
symmetry to reduce the squared-error to a sum of the form
$E_{\lambda_{1},\dots,\lambda_{d}}\sum_{i}g(\lambda_{i})$. The main difference
is that in this case the eigenvalues of $G$ are no longer independent, but
they are at least still exchangeable, so we can bring out the factor of
$|\beta_{0}|^{2}$ in front of the sum like above. And since $G$ now has a
Wishart distribution, we can use the Marchenko-Pastur theorem (Bai &
Silverstein, 2010) instead of the Law of large numbers to reduce the sum to
the indicated integral form.

###  A.7 Estimation of minima and curvatures

To estimate the minimum and curvature of an error curve $Err(\alpha)$, we
evaluate $Err(\alpha_{i}),i=1,\dots,n$, where where $n=500$ and $\alpha_{i}$
are equally logarithmically spaced between $10^{-3}$ and $10^{5}$. The minimum
is simply estimated as $\min_{i}Err(\alpha_{i})$. To estimate the curvature at
the minimum, we took the closest few points to the minimum and fit a quadratic
function. That is, if $i_{0}=argmin_{i}Err(\alpha_{i})$ then we fit a two-
parameter linear model of the form

$Err(\alpha_{i})-Err(\alpha_{i_{0}})\sim
a(\alpha_{i}-\alpha_{i_{0}})+b(\alpha_{i}-\alpha_{i_{0}})^{2},i_{0}-5\leq
i\leq i_{0}+5$ (68)

with the estimated Hessian being $2b$.

###  A.8 Relation between $\alpha$ and $C$ in Theorem 2

To express the relation between $\alpha$ and $C$, we consider two cases,
corresponding to whether or not the constraint set contains the global
minimizer of the objective $L\mapsto Tr(LL^{T})$. The first case is when
$C\geq\|I_{d}\|_{p}=d^{1/p}$. In this case, the optimum is evidently
$L=0_{d\times N}$, corresponding to $\alpha=\infty$.

In the case where $C<d^{1/p}$, $0_{d\times N}$ does not lie in the constraint
set. Therefore, the optimum $L$ lies on the boundary of the constraint set,
i.e. $\|LX-I\|_{p}=C$. By plugging in the functional forms from Theorem 2, we
obtain the following relations:

$\displaystyle\sum_{\sigma_{i}<\alpha}1-{\frac{\sigma_{i}}{\alpha}}$
$\displaystyle=$ $\displaystyle C\hskip 28.45274pt(p=1)$ (69)
$\displaystyle\sqrt{\sum_{i}{\frac{\alpha^{2}}{(\sigma_{i}+\alpha)^{2}}}}$
$\displaystyle=$ $\displaystyle C\hskip 28.45274pt(p=2)$ (70)
$\displaystyle\alpha/(1+\alpha)$ $\displaystyle=$ $\displaystyle C\hskip
28.45274pt(p=\infty)$ (71)

where $\sigma_{i}$ are the eigenvalues of $X^{T}X$.

###  A.9 Effect of number of $\alpha$ values

As pointed out in the main text, the performance of the cross-validated models
can depend on the number $n$ of $\alpha$ values used for cross-validation. To
do so, we follow the methodology used in 3.2.1, with the only difference being
we use a smaller hyperparameter range
$\sigma\in[.5,2,3.5],\lambda=.5,\rho\in[0,.5,.9]$, and also vary the total
number $n$ of $\alpha$ values used in the cross validation
$n\in[9,15,20,30,50]$. We keep the limits of the range of $\alpha$ values as
before, and also maintain the equal logarithmic spacing.

We show the average error and win probability in Figure 7 and Figure 8
respectively. As per the discussion in Sections 2.2.3 and 5, we see that the
Ridge often benefits drastically from increased number of $\alpha$s, and can
overtake the Nuclear for large number of $\alpha$ in cases when the Nuclear
performs better for small number of $\alpha$.

Figure 7: Average test error for varied number of regularization strength
values $\alpha$ used in cross-validation. Each point represents an aggregate
over 100 datasets.  Figure 8: Same as Figure 7, except showing the probability
that each model attains the lowest test error on a given dataset. Each point
represents an aggregate over 100 datasets.

###  A.10 Sparsely structured data and comparison to Lasso

We consider a variant of the setup in Section 3.2.1, in which the data is
constructed to have sparse structure. In this case, when generating the ground
truth coefficient vector $\beta_{0}$, we select a set of indices
$I\subset\\{1,\dots,10\\}$, $|I|=3$ at random, and generate $\beta_{0}$ as

$\displaystyle(\beta_{0})_{i}$ $\displaystyle=$ $\displaystyle N(0,1),i\in I$
(73) $\displaystyle(\beta_{0})_{i}$ $\displaystyle=$ $\displaystyle
N(0,1)/10,i\not\in I;$ (74)

We also use the smaller hyperparameter ranges
$\rho=0,\lambda=.5,\sigma\in[.5,1,1.5,2,2.5,3]$. Otherwise we follow the
methodology of Section 3.2.1.

We also include Lasso in the set of considered models, since it is designed to
deal with sparse coefficient vectors. Note, however, that Lasso is not
manifestly a Linear Estimator 101010We do not, however, have a proof that it
cannot be expressed in the requisite form of Equation 1 in the sense defined
in Section 2.1 .

We show the average error and win probability in Figure 9 and Figure 10
respectively.

Figure 9: Average test error, with 95 percent confidence intervals, on
Gaussian data with sparse ground truth coefficient vector. Each bar is an
aggregate over 100 datasets. Figure 10: Same as Figure 9, except showing the
probability that each model attains the lowest test error on a given dataset.
Each bar is an agggregate over 100 datasets.

###  A.11 Real data experiments

We evaluate the models on real (i.e., non-synthetic) data. We consider the
well-known Diabetes dataset (N=442, d=10) and California housing dataset
(N=20640,d=8), both available from the `sklearn.datasets` library.

To analyze the performance of the models on each dataset, we create random
train-test splits in which the size of the training set is always set to 300,
and the test set comprises the remaining observations. Each model is fit on
the training set (including the regularization strength $\alpha$, using the
same cross-validation procedure as in Section 3.2.1), and the mean-square-
error is evaluated on the test set. We use the same set of 9 $\alpha$ values
as in Section 3.2.1 for all models. We constructed a total of 200 splits in
this way , and evaluated each model on each split as before (so that we can
evaluate each model’s probability of winning on a given split, as well as the
average error over splits).

We show the average error and win probability in Figure 11 and Figure 12
respectively. We see a similar pattern as in Figure 5, with the Nuclear having
a clear advantage over the other models in terms of win probability. The
Nuclear also attains the lowest average MSE in the diabetes data, and is
essentially tied with Ridge for lowest average MSE in the housing data.

Figure 11: Average test error, with 95 percent confidence intervals, on a
random train-test split, for the diabetes and California housing datasets.
Each bar is an aggregate over 200 splits.  Figure 12: Same as Figure 11,
except showing the probability that each model attains the lowest test error
on a given split. Each bar is an aggregate over 200 splits.