Flat Minima in Linear Estimation and an Extended Gauss Markov Theorem
Authors: Simon Segert
Venue: ICLR 2024 poster (2024)
Keywords: regularization
regression
optimization
Primary area: learning theory

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.

Published as a conference paper at ICLR 2024
FLAT MINIMA IN LINEAR ESTIMATION AND AN EX- TENDED GAUSS MARKOV THEOREM
Simon Segert Princeton Neuroscience Institute Princcton University
Eneegert@prinnnetoo edd
ABSTRACT
We consider the problem of linear estimation, and establish an extension the Gauss-Markov theorem, in which the bias operator is allowed 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.
INTRODUCTION
Linear models are among the most used of all machine-learning models in applications to science and engineering. In addition this practical interest, they 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 (heory(Belkin al., 2019p. This has begun to be reconciled by careful analysis of well-chosen linear "model systems" such as high- sigh-dimensionall regression (Hastic et 2019), kernel ridge regression (Canatar et al.] 2021: Jacot al., 2020). or linear neural networks (Saxe et al. 2013), which can reproduce many qualitative properties of learn- ing 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 key ingredient for achieving high performance of both linear and non-linear models. For example, neural networks are often
trained using L2 regularization (Goodfellow al., 2016) dropout (Srivastava al. 2014 By contrast, most of the above-mentioned thcoretical 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 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 (Kamalak lakara et al. 2022 Yang
al. 2020) Spectral norms (Johansson et al. 2022 Yoshida & Miyato 2017). The Nuclear norm (or convexified rank) also plays key role in other kinds of bighddimensional non- non-linear learning problems (Hu et al. (02T), most notably matrix completion, and is related to dropout regularization
in neural networks (Mianjy & Arora 2019).
Thus with the general motivation of gaining detailed understanding of different kinds of regular- izations in a tractable yet reviiuuslyvvalddated and informative setting, 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 Gauss-Man Markov-likk theorem for this setting, which characterizes the optimal estimator at 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.

---

Published as a conference paper at ICLR 2024
Next, 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 always attain the lowest enror of any estimator, the minima
pattern: while the Ridge estimator can be sharp (as a function of regularization strength). The Nuclear estimator, by contrast, can often attain test error nearly low as the Ridge case, but with significantly flatter minima. This has implications for practical situations when the regularization strength is selected using noisy search procedure such as cross validation, as could happen that such procedure finds better solution in the latter case, even if the actual global optimum is better in the former case. We then perform series of simulation studies with cross-validese validaaai,,,,nn, 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-MARKON THEOREM FOR SCHATTEN NORMS
Throughout, we will be concerned with the classic regression problem. Denote the training data by X, which is an N x d matrix, and the training targets by Y. which is N x vector. Y is assumed to have distribution Y = X + c, where c are independent samples from a fixed noise distribution. We will assume that the noise distribution has mean zero and variance 0² 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
3 LY
(1)
Where L is an d x matrix that depends (possibly non-linearly) on the training observations X. Given a new datapoint the predicted value is just ÿ = (Ztest: 3). One simple observation is that for any linear estimator the bias E, - 3u can be expressed as a linear function of 3p:
E, B 3o (LX 1)
(2)
We thus define the bias operator B == LX I. Note that the actual bias value generally depends on the unknown quantity fo. while the variance Varz( L.I does not. The classic Gauss- Markov theorem tells us that if want to impose exactly zero biad" 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 Sg. As a proxy, we propose to control the size of B instead, using some choice of matrix norm. what follows, we will only consider Schatten norms, 2 1, although in principle other matrix norms could be used. This
|| = (E, o, (M)))
motivates the following:
dd the p-Bias constrained Linear estimator (X) is defined as the minimal variance estimator which has p-bias of at most C. That is,
Definition 1. Let p 2 1, and let C 220 O. For any matrix X E N
Lx(X)= brrgninne***,jixxxil.ce
where I1p is Schatten norm. Y a vector of regression targets, the estimated coefficient vector is = (X)Y.
We now state our main theorem:
Theorem 2. (Extended Gauss-Markov) Let X Nxd be matrix of observations, Assume that G= XT is invertible, and let G = Udiag(o) U be the diagonalized form. For p 2 1, the p-Bias constrained estimator with bound C' can be expressed in the form Lp (X) = G where & is coametrrc,simultaneeee cases of p we further have:
diagonalizable with G, and satisfies G c For the following special
Note the subtic point that this is the only case in which we can exactly control the bias without knowing the value of Bu
We reserve the notation || ||p for a Schatten norm of matrix, and use Ip for Euclidean vector norms
i.c. G
G non-negative definite
2

---

Published as a conference paper at ICLR 2024
p 1 (Nuclear norm): G ['dio a))
p 2 (Frobenius norm): G G ala
p oc (Spectral norm): G (1 a)G
where a 2 determined by C, and max denotes the elementwise maximum, i.e. max a); = mar(a;; a)
Note that the p 2 case is just Ridge regression, while the p case is just scalar shrinkage of Bor. towards 0, as per Stein 1962). We refer to the p 1 case Nuclear Regression and likewise Spectral Regression for the 2o case. The exact relation hetween C' and a is given in A.8
We also note that theorem does not perfectly generalize the classical Gauss-Markor theorem. Var(ß) where 8 is any unbiased linear estimator (Wooldride 2015). However, taking C' = 0 in our Theorem we would conclude only that Tr(Van a (BooL))) S Tr(Var(3) which is strictly weaker.
Indeed, the classical theorem is usually stated in the form Var(Sor.s)
2.5
ridge
2.5
ridge nuctear
2.5
rige nuclear
nucinar
1.5
0.0
0.0
0.0
0.0
U.5
1.0
4.5
2.0
0.0
0.2
0.4
0.6 U.M 1.0
Figure 1: Illustration of Theorem 2 Each point the curve corresponds to regression fit with certain value of the matrix X Viag(VI, v10). We show the norm of the bias matrix B LX - on the axis, and the sample variance Tr(L" L) on the y axis. The green curve always lies below the other two the left plot, and analogously for the other two plots.
2.2 FORMULAS FOR TEST ERROR
We next analytically characterize the generalization performance the estimators derived in the previous section, under certain specific assumptions on the distribution of X. We will consider two diffcrent 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, an appropriate thermodynamic limit.
2.2.1 SPHERICAL GAUSSIAN ENSEMBLE
Consider predictor matrix X Nxd, where each entry is an independent centered Gaussian of variance 1/N. We consider the thermodynamic limit in which N > oo and d/N + A < 1. The training targets are distributed as Y, (X,, 30) + oee, where E; N(0, 1), and the magnitude of the ground-trus trutt coefficient vector fo satisfies 1Bo12 + 32, for some fixed B > 0. The testing error for given dataset X, Y is defined as
MSE == E ((x, Bo) - (x, B))
(4)
where 3 is the coefficient vector estimated from X and Y. (Note that we do not include exogenous noise the testing data; if we did, would just contribute an additive offset of g²). We observe that the expression can be simplified to 1.3g - /12/N, however we write it in the more verbose form above to make clear that is obtained 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 show (cf. A.6):
Proposition 2.1. Consider the p-bias constrained estimator with regularization strength
3

---

Published as a conference paper at ICLR 2024
In the above limit, the average test evvor Erry(a) Ex,yyyMEEE.. is given by I
Err(a)= (am(a)=>J (81 fa(x) fa(o))/Mroddd)
0² x
(5)
where PM is the Marchenko-Part Pastur density with concentration parameter X, and fa is defined casewise as follows:
[a(x) maz(s, p), fa(2) a. p
(7)
fa(x) = = oc
(8)
In practice, we use this formula and evaluate the integral using numerical quadrature. However, is interesting note that the integral 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 HMP P that easily follows from the above integral (Bai & Silverstein 2010 In the spectral case, a formula using simple calculus:
derive
Proposition 2.2. For the Spectral estimator in the spherical Gaussian setup, the test enror in the thermodynamic limit equals
A3 + 1-
Err
(9)
a1+
To state the result for the Nuclear case, is necessary to introduce the Appel hypergeometric function Fi (Applll 1925p, 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 functioo Fi:
a
Erri(a)= CD + Va-a - (1) 2.2 x
X4 X4-
1_ SaSA
Ao
-2022 + X32
a a d+
a²
where X are the limits of the Marchenko- Pasttur support: = (1+VX)2 CDFx is the cumulative distribution function of pMr and ca((a) are certain rational functions of
See A.5 for proof, and precise definition F1-
2.2.2 DIAGONAL ENSEMBLE
While the spherical Gaussian case is an instructive starting point, the assumptions are clearly some- what limiting. In particular, it implies 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 essentially arbitrary spectral density. On the one hand, this new ensemble has its own limiting assumptions that the Gram matrix X is almost surely diagonal, longer independent. But the other, having
and the observations (i.e. rows of X) other random matrix model allows us to assess whether qualitative observations about the spherical Gaussian case generalize other settings.
We now describe the random matrix ensemble. We first specify some distribution v supported on [0, 1], which will be the spectral density. We also specify some distribution Vnoise supported on bounded interval of Rt, which will act as multiplicative random noise on the spectrum; we require that Er.w.... F To generate training and testing data, we first sample X; u, 1,...,d independently, and 8; Vnoise- We then set
X, = Xidigg( (pst))); Xtest X2diag( As)e)
(10)
where X, and X2 are drawn independently and uniformly from the Stiefel manifold O(n, d) of N x orthogonal matrices. That is X X; = la- We observe that, as noted above: 1.) the Gram matrix G X X,r diagonal almost surely, and 2.) the limiting eigenvalue density of G is given by
4

---

Published as a conference paper at ICLR 2024
The regression targets are then constructed as before: Yo. Xe + c, Yeest Sve The test error is defined as
MSE = N |X Xtest/0l
(11)
Note that this is slightly different from the expression for MSE in the spherical Gaussian case (Equa- tion The basic reason 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 time and average, here we sample batch of N (dependent) test observations a time, and compute the aver- age eπrоr 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 (r) = X MMM. SE can be expressed as
Errpla) S Px(1
x² La Ja(
2
(12)
(x)
where Ja has the same form as in Proposition 2.7
See A.6 for proof. Note the similarity to the result of Proposition The basic reason for this is that both ensembles have a similar rotational invariance, which allows the MSE be expressed as sum over the spectrum of Xer. 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 dv(x)/dz x x 10. as many
real covariance matrices have approximately this structure (Stringer et al. 2019; Qin & Colwell, 2018 Liu et al. [197] Ruderman & Bialek, 1994p. In this case, the expressions for the test error can again be simplified to analytic expressions that do not involve integrals, see Section A.9 for details.
2.2.3 ON THE OPTIMAL REGULARIZATION FORM
It not hard see that the derivations in the above sections go through for any estimator of the form 3 Y-T, where the eigenvalues of G are related those of G through some fixed point-wise transform (G) (G), and & 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 zero it is easy to see that in both of the above random matrix models, the optimal form of corresponds to Ridge regression with regularization strength o² /32.
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 3 and a, which typically not known. In practice optimal ridge parameter us..g.e.g.. cross-validation, and hope that our estimate is close to the oracle- optimal value d²/ In the finite-sample regime, there will be some noise in the cross- cross-validated errors, so this essentially amounts to trying to minimize Err(a) given only access to 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 small number of values of a. very rough rule of thumb that if we sample n values ag all within some small distance (r66²
would estimate the
8 to the minimum, then min, Err(os) R f + [n+1][++2)' where minimum, and x² is the Hessian of Err(a) at the minimum (see 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.
min, Err(x) is the true
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 and compared the average test error with the the- oretically predicted values. For each such dataset, we used = 100 training observations, and obtained the empirical test error by averaging over 100 datasets. We did this for both the spherical
this be justified using results from Calculus of Variations
5

---

Published as a conference paper at ICLR 2024
Gaussian case, well as the diagonal case with power law spectral density dv/dz x x>> 1nc (o. In Figure 2 we show several plots with different values of the other hyperparameters, in all cases we 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 often noticeably flatter than those of the Ridge, while also being only moderately higher.
Am
ridge
riöge
spectral
spectral
0.36
0.32
0.08
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.Not that the minima for the Ridge curve are notably steeper than for the other two.
EXPERIMENTS
3.1 LOSS BASIN GEOMETRY COMPARISON
Here, we aim systematize the observation that the Nuclear estimator have slightly higher, but much flatter, minima than the Ridge estimator. To do so, we simply fix the hyperparameters a, A, and consider the theoretical test errо as a function a. We record both the minimal value over a m = Eir(a) as well as the curvature at the minimum: == ((PErr(())/ee la=amsn. For all cases, we fix 3 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 For the nuclear case, we can plainly see that the minima of the test loss are nearly as deep for the the Ridge (within 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 less favorable depth tradeoff compared to the Nuclear. Moreover, we see a pronounced effect of ar, 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 we see that these observations largely carry over to the case of power-law spectral density, indicating that they model. See for more technical details on the calculation of these values.
simply incidental property of the spherical Gaussian
nuclear
minicuvatuue min, relative to Ridge
spectral
5-5.2
A-04
A-34
1:04
A-12
2:0.4
A-24
2:04
P-C.5
12.7/91
21.25535
12:2/3255
1248.3
18024.3
b.405'
2.828
0.6-59
03333
5.4020
73646.
33.2165
51255
03:30
0153
0.145
092.5
2,8-22
$3/472
7,9989
108/971
Figure 3: Quantification of loss basin geometry for the spherical Gaussian model. Each entry 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.

---

Published as a conference paper at ICLR 2024
nuclear.y=..3
spectral.ym0.5
P-0.25
12.3/70.7
19.3/-6.7
12.3:76.7
12.3/77..7
9-0.25
26.0/422.
26,0/62
26.0/42.
26.0/662.8
oe0 75
3.1/9.5
3.1-8.5
3.1/6.5
3.1/8.5
dw0.75
22.2/50..4
22.2/:004
22.2/-0.4
22.2/-0..4
n
spectral.y
Aw0.2
Aw0.4
Awo
Awd.8
Aw0.2
Aw0.4
Awd
A=D
0-0.5
11:5/44.6
11:5/492.0
11:5/42.0
11:5/42.0
0-0.5
11:5/434
11:5/-144
11.5/444
11. 11.4/43.4
0.9/9.2
0.9/0.2
0.9/0.2
0.9/0.2
8.6/49.5
5.6/4855
0.61-9.5
0.8/40.5
nuclear.y=1..5
spectral.y=1.5
Λ-0.6
A-0.8
A=0.2
Λw0.6
10.4/:00.5
10.4/80.5
20.4/800.5
10.4/55..5
de0.25
3.3/42.3
3.5/422
3.5/422.
3.5/42.3
10.1/:5077
0-0.75
2.5/1.1
255-111
2.5/1.11
2.5/1.1
9-0.75
6.27-228
6.2:-128
6.21-428
6.2/212.8
Am0.2
Aw0.4
Awo.6
A-0.0
Am0.2
Aw0.4
Awo.6
A-0.0
9.7/-03.9
0.7:03.9
0.7/-03.9
1.9/35.7
1.9/-35.7
1:9/-3577
1.9/35.7
dHw0. 5
8.7/225.2
8.7/25..
6.71-25.2
8.7/225.2
de0.3
3.7/-515
5.7/31
3.7/-31.5
3.7/31.5
p=1
0.7//0.7
0.716.7
0.7/0.77
0.7/0.7
g=1
3.5/31.6
2.9/11:6
3.5/-31.6
2:5/-116
Figure 4: Same as Figure B except for the Diagonal matrix model with power law spectral density.
3.2 SIMULATION STUDIES
The purpose of this section 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. GAUSSIAN PREDICTORS
We generated training data by sampling X, ~ C(p)), where (1 pl1'. We then generated ground truth coefficient vector by sampling Bo N(0, Ia), and training targets by N((X,, Bo), o²). The testing data were generated similarly. using the same $o. In all cases, we used N 100 training examples, and 5000 testing examples. For each combination of hyperparameters d, a. 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 best-performing a on the training set. We then refit the model the entire training set, and evaluated its performance on the test set. The set of allowable a values considered the cross-validatio was the same for all models, and consisted of 9 equally logarithmically spaced values spanning 10-1 to 106.
In Figure 5 we plot the best-performing model for each combination of hyperparameters, for two different definitions of "best". In the top row, 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 test error of M EmE: The top row shows We see that in terms of average error, the Ridge and Nuclear are essentially matched, with Spectral being distant 3rd. The Nuclear however is overwhelmingly likely to be the winner on given dataset.
MS
while the bottom row shows M odee gargminm.MM..
hy).
3.2.2 REGRESSION 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 kemel. To be more specific, the process for generating data was follows. We first pick integers d, dras which

---

Published as a conference paper at ICLR 2024
pe50
d-00
223
PROK
8
54025
5
1:20, pob best
6-42 pred test
0-60, pab best
(-31 prat
F213
2922
3-30
2+L3
24
2+
20
2222
0425
$-35
0413
2V
2425
80
15
8 3
7 3
2-33
11
3
3
2403
6455
Figure 5: Best performing models with cross- volssvvaldattinn Gaussian predictors. Top: best average emor. Bottom: highest probability of winning. Each cell an aggregate over 100 datasets.
define the dimension of the observed space, and of the kernel feature map, respectively. We sampled training matrix X of size N x d, in which each entry has variance 1/d,of. We then sampled non-linear function f I + IR in a manner described below, and generated training targets as y: = f(X,)++++) to,, where X, the ith row of (i.e. ith observation) and e; is unit Gaussian. Testing data was generated similarly, using the very same function (and setting analyses). Since the relation between X and Y non-linear, we do not directly fit the models on X, but rather first fit Random Fourier features (RFF) model (Rahimi & Recht 2008) drof features We then fit the linear models on F(X), ) Y, where RFF(X) is the N dybf matrix obtained by applying the Random Fourier features mapping. To estimate test error, first transform test using the same RFF mapping, and consider [Xtes())))) Yeese Finally, the non-linear function f was defined as follows:
0 as in the previous
X, using
drss
f(x)= = Σ (2 wx))/k²
(13)
where are sampled uniformly from the unit sphere 5 C Rdd We fixed d = 10 and 100 all cases. For each value of drss and o, 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 As in Figure [ 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 a.
1.2
2,0
1.0
0.75
nucCee
0.70
0.65
0.7
0.5
0.80
1.0
0.55
0.30
0.5
200
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 average Ridge test error. Right shows probability of "winning" on a given dataset. Each dot represents an aggregate over 100 datasets.
"The RFF fit done using the class sk eeeearnn.. kernel 1_aprne_ppppxxiiii RHM'Smp er.
8

---

Published as a conference paper at ICLR 2024
4 RELATED WORK
Several other sorts of generalized Gauss-Markov theorems have appeared in the literature, although they generalize along different axes than our formulation, such as allowing non-spherical errors (Kourouklis & Paige 1981). predictor collinearity (Lewis & Odell I966p, random effects (Shaffer 1991]. or non-lineat estimators (Hansen 2022). Similarly,the Nuclear and Spectral estimators bear a very strong resemblance to various notions from classical statistics, most notably Principal Com- ponents Regression (PCR) (Kendall 1957), Stein Shrinkage (Stein) (162)) respectively. Finally, the work of Hocking et al.] (176) is particularly relevant and worthy of further discussion here-we were-wegge give more detailed comparison in A.2
Lately, there there has also been renewed interest 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 al. 2023) and early stopping (Ali et al. 2019). This pervasiveness of the ridge underscores our claim that 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 (201 I); Samet al. (2013) Lasso. Theoretical properties of nuclear norm have mainly been studied in the context of matrix completion, in which it is often used as surrogate for matrix rank (Candes & Recht 2012 Koltchinski al. 2011). but see Hu et al. (202T) 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 has to do with regulariz- ing the outputs of a regression problem with multiple dependent variables. There is also some work There
on theoretical properties of the Principal Components Regression estimator (Xu Hsu 2019;; opti- mizing the number of PCR components has similar flavor to optimizing in the Nuclear estimator.
Finally, flat minima are also popular as an explanatory tool in the context of generalization of neural networks (Baldassi et al.) | Mulayoff & Michaeli 2020 Hochreiter & Schmidhuber] 1997). in which the randomness in the minimization procedure is induced by stochastic gradient descent.
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 a val- ues 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 are used. However, increasing the number of values has other potential costs such as risk of overfitting or increased computational cost. But more fun- damentally, 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 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, hope that this framework will provide test case understand less common regularizations such as || l| and || Il with the goal of transferring the insights gained thereby to the study of high-dimensional non-linear models.

---

Published as a conference paper at ICLR 2024
ACKNOWLEDGMENTS
SS is supported by a T32 training grant in Computational Neuroscience (T32MH0652144)
REFERENCES
Ali, JZ Kolter, and RJ Tibshirani. A continuous-tima view of early stopping for least squares. AISTATS, 2019.
Appell. Sur les fonctions hypergéométriques de plusieurs variables. Mémoir Sci. Math. Paris: Gauthier-Villan 1925.
Bai and Silverstein. Spectral Analysis of Large Dimensional Random Matrices, Springer, 2010.
WN Bailey. On the reducibility of appell's function f4. Quart. J. Math, 1934.
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.
M Bayati and A Montanari. The dynamics of message passing on dense graphs, with applications to compressed sensing. IEEE Transcations on Information Theory, 2011.
M Belkin, D Hsu, Ma, and S Mandal. Reconciling modern machine learning practice and the bias-variance trade-off. Proc Nat Aca Sei, 2019.
M Belkin, D Hsu, and J Xu. Two models of double descent for weak features. preprint arXiv: 1903.0751, 2020.
F Bunea, Y She, and M Wegkamp. Optimal selection of reduced rank estimators of high-dimensional matrices. Ann. Starist., 2011.
A Canatar, B Bordelon, and C Pehlevan. Spectral bias and task-model alignment explain general- ization kernel regression and infinitely wide neural networks. Nature commumications, 2021.
E Candès and B Recht. Exact matrix completion via convex optimization. Communications of the ACM, 2012.
LH Dicker. Ridge regression and asymptotic minimax estimation over spheres of growing dimen- sion. International Statistical Institute and the Bernoulli Society for Mathematical Statistics and Probability, 2016.
Dobriban and S Wager. High-dimensional asymptotics of prediction: Ridge regression and clas- sification. Annals of statistics, 2018.
Ky Fan. Maximum properties and inequalities for the eigenvalues of completely continuous opera- tors. Proc. Nat. Acad. Sci. U.S.A., 1951.
I Goodfellow, Y Bengio, and A Courville. Deep learning. MIT press, 2016.
BE Hansen. A modern gauss- markov theorem. Econometrica, 2022.
T Hastie, A Montanari, S Rosset, and RJ Tibshirani. Surprises in high-dimensional ridgeless least squares interpolation. srxiv.abs/1003.08555 2019.
S Hochreiter and Schmidhuber. Flat minima. Neural Computation, 1997.
R. R. Hocking, F. M. Speed, and M. Lynn. A class of biased estimators in linear regression. Technometrics, 1976.
H Hu and YM Lu. Sharp asymptotics of kemne ridge regression beyond the linear regime. 27Xiv.2205.06798, 2022.
Z. Hu, F Nie, R Wang, and X Li. Low rank regularization: A review. Neural Networks, 2021.
10

---

Published as a conference paper at ICLR 2024
AJ Izenman. Modern multivariate statistical techniques: Regres- sion, classification, and manifold learning. springer, 2008
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, Engsner, Strannegard, and P Mostad. Exact spectral norm regularization for neural networks. arxiv:2206. 1at.e220...1381.... 13581.pdf. 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.
MG Kendall. A Course in Multivariate Analysis. London: Charles Griffin, 1957.
V Koltchinskii, K Lounici, and AB Tsybakov. Nuclear-norm penalization and optimal rates for noisy low-rank matrix completion. Ann. Statist, 2011.
S Kourouklis and CC Paige. A constrained least squares approach to the general gauss- markov linear model. Journal of American Statistical Assoc., 1981.
TO Lewis and PL Odell. A generalization of the gauss-markov theorem. Journal of American statistical assoc., 1966.
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, P Cizeau, M Meyer, CK Peng, and HE Stanley. Correlations in economic time series. Physica A: Statistical Mechanics and its Applications, 1997.
S Mei and A Montanari. The generalization error of random features regression: Precise asymptotics and double descent curve. arxiv.org/pdff119 http://///////....ooo http:////00000055...... 2020.
P Mianjy and R Arora. On dropout and nuclear regularization. International conference machine Learning, 2019.
R Mulayoff and T Michaeli. Unique properties of flat minima in deep networks. Proceedings of the 37th International Conference on Machine Learning. 2020.
Nakkiran, Venkat, SM Kakade, and Ma. Optimal regularization can mitigate double descent. In International Conference on Learning Representations, 2021.
Qin and LJ Colwell. Power law tails in phylogenetic systems. Proc. Nat. Aca. Sci, 2018.
Rahimi and Recht. Weighted sums of random kitchen sinks: Replacing minimization with randomization in learning. In NeurlPS, 2008.
DL Ruderman and W Bialek. Statistics of natural images: Scaling in the woods. Physical review letters, 1994.
Samet, C Thrampoulidis, and B Hassibi. The squared-error of generalized lasso: A precise analysis. Allerton Conference on Communication. Control, and Computing. 2013.
AM Saxe, JL McClelland, and Ganguli. Exact solutions to the nonlinear dynamics of learning in deep linear neural networks. arXiv:1312.6120) 2013.
JP Shaffer. The gauss markov theorem and random regressors. The American Staristician, 1991.
Srivastava, Hinton, Krizhevsky, Sutskever, and R Salakhutdinov. Dropout: simple way to prevent neural networks from overfitting. The Journal of Machine Learning Research, 2014.
J Stein. Multiple regression. Contributions to Probability and Statistics. Essays in Honor of Harold Hotelling, pp. 424-443, 1962.
Stringer, M Pachitariu, N Steinmetz, Carandini, and KD Harris. High-dimensional geometry of population responses in visual cortex. Nature, 2019.

---

Published as a conference paper at ICLR 2024
S Wager, S Wang, and P Liang. Dropout training as adaptive regularization. arXiv: 2013.
JM Wooldride. Introductory Econometrics: A Modern Approach. Cengage Learning, 2015.
J Xu and D Hsu. On the number of variables to use in principal component regression. Neurips, 2019.
Yang, M Tang, W Wen, Yan, Hu, Li, H Li, and Chen. Learning low-rank deep neu- ral networks via singular vector orthogonality regularization and singular value sparsification. arXiv:2004.09031, 2020.
Yoshida and T Miyato. Spectral norm regularization for improving the generalizability of deep learning. arXiv: /705. 720./7055.1490 10941, 2017.
12


Reviews:
[{"review_id": "dohISJ0ASE", "reviewer": "Reviewer_cHcu", "summary": "The Gauss-Markov theorem says that the best unbiased linear estimator is the pseudoinverse of the data matrix (in the paper denoted by $X$), which can be obtained by minimizing the Frobenius norm of the estimator $L$ subject to the constraint that $L$ is the left inverse of the data matrix, $LX=I$. The paper generalizes this formulation by relaxing this constraint to $||LX-I||_p\\le C$ where $||\\cdot||_p$ is the Schatten $p$-norm (vector $p$-norm on the singular values), thus allowing some bias $C$. A special case is the ridge/Tichonov regression, obtained for $p=2$. The parameter $C$ (or a monotonically related parameter $\\alpha$) is determined by validation, by minimizing the test error. The advantage is that for $p\\neq2$, the minimum of the test error over $C$ may be flatter than for $p=2$. \n\nThe authors derive a closed form solution for the optimal estimator $L$ for $p=1$ (nuclear norm) and $p=\\infty$ (spectral norm). Next, they derive explicit form (as integrals) of for the test error in thermodynamical limit for two special distributions of the data: spherical Gaussians (elements of $X$ and noise in right-hand sides are i.i.d. normal) and diagonal data (when the Gramian $X^TX$ is diagonal). The integrals are solved in closed form for $p=1$ and $p=\\infty$. This theoretical formula is shown to agree with test error on synthetically generated data.\n\nThe cases $p=1,2,\\infty$ are compared in a simulated experiment in which $\\alpha$ with smallest test error is estimated by x-validation, where best $\\alpha$ is selected by \"grid search\" from 9 log-spaced values. This showed that the nuclear-norm regression is comparable to or better than (depending on methodology) ridge regression. A similar result is obtained for the similar experiments for nonlinear regressors (random Fourier features).", "questions": "It would be helpful to vary the number of values of $\\alpha$ in the grid method (currently this value is 9) in the experiments, i.e., to make it a hyperparameter. Pls see my remark on this in \"weaknesses\".\n\nIt would be helpful in Figure 5 to report not only winners but also MSE for different models (as in Figure 6 left) - because a winner can win by only a small margin.\n\nWhy are any experiments on real data not included? Is there a theoretical obstacle? Or, perhaps, you believe that the results of such experiments would not be informative enough? Please comment (also in the paper, if accepted).", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The observation that a wider minima of $\\alpha$ can be achieved at the cost of a little bias is interesting.\nThe theoretical results (test errors) are non-trivial to derive.\nThe text is clear enough, though clarity could be improved by more effort.", "weaknesses": "I cannot assess novelty reliably because my expertise is mainly in optimisation rather than estimation (however, I understand all parts of the main paper well). In fact, rather than extending Gauss-Markov theorem, the paper generalizes ridge regression (please, consider changing the title accordingly). It is well-known that ridge regression has non-zero bias but a smaller variance than pseudoinverse.\nSo the main novelty seems to be that of flatter minima of test error, rather than the generalization of ridge regression.\n\nThough the minima for the nuclear-norm regression are indeed flatter than for the ridge regression, the difference is sometimes only minor, as seen in Figure 2. It is true that the experiments show that in estimating $\\alpha$ by x-correlation, the nuclear norm most often wins. However, this might be due to the experimental methodology. E.g., if there were more than 9 values of $\\alpha$, the deeper minima of the ridge regression might have been hit much more often.\n\nA major weakness, in my opinion, is that the experiments are done only on synthetic data. The applications of linear regression are abundant, so it should be possible to find many suitable real datasets for this.\n\nMinor/fixable issues:\n\n1st formula in section 2.1: symbol $L(X)$ is used here but then never more. Change to $L$.\n\nThm 2: Letter $\\Sigma$ is usually used to denote the diagonal matrix with singular values. For vector of singular values, better use $\\sigma$ or $s$.\n\nThere are many small mistakes/typos in the text. E.g., references to figures in sections 2.2.4 and 3 refer to non-existent figures (e.g., figure 2.2.4 in section 2.2.4).\n\nIn the 5th line of section 2.2.1, can the formula for MSE be simplified (i.e., calculate the mean value in closed form)? It is confusing that the MSE has quite different form in sections 2.2.1 and 2.2.2 (this lets the reader wonder if this difference is substantial or just due to little care for text clarity). This might deserve a comment.\n\nMost displayed equations are unnumbered, which is not friendly for reviewers (given that lines are not numbered in the ICLR style).\n\nThe asterisk symbol in Proposition 2.1 and in the first displayed formula in section 2.2.2 has not been defined.\n\nThe integrals for $Err(\\alpha)$ in sections 2.2.1 and 2.2.2 are almost the same, up to the integrating measures. I wonder if they are correct or there are typos in them..?\n\nPOST REBUTTAL: The authors have clarified my objections, to certain extent. I am therefore raising my evaluation. However, please note that I cannot assess novelty reliably.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "TDtw25ZBrf", "reviewer": "Reviewer_ACbc", "summary": "The paper considers a variant of linear regression with constraints placed on a \"bias operator\". Under this framework, the paper discusses an extension of the Gaussian-Markov theorem, showing empirical and theoretical evidence for its main result, Theorem 2. Later discussions in the paper surround \"flatness\" and \"deepness\" of various losses considered within the paper.", "questions": "Several questions were listed in the above \"weaknesses\" section.", "limitations": "", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The paper tackles an interesting class of linear regression models and delivers a thorough presentation of various aspects of the problem, from the problem definition, main theorem statement, to several case studies, all of which help to paint the overall picture of the problem. The constrained setup considered in the paper is also interesting and intuitive. Considering that linear models are a core concept of machine learning, the paper is of sufficient interest to ICLR.", "weaknesses": "There are several dimensions of weaknesses presented in the paper:\n\na. Clarity and overall quality of presentation. The paper does not appear to be carefully edited and revised, with multiple typographical errors in the first paragraph of the introduction alone (examples: \"somewht\" in line 3, lack of period at end of sentence in line 4, reverted quotation marks on line 5, etc.). The graphs, equations, and tables in latter parts of the paper can also benefit from detailed revisions. These issues surrounding clarity and presentation are not constrained to the first paragraph and can be found throughout the paper and also the supplementary material.\n\nb. Discussion of main theorem. While the problem setting itself is interesting, the discussions of the main theorem (Theorem 2) leave an impression that it can be further discussed. For example, what are the values of alpha? Although it is shown in the supplementary that alpha is a consequence of solving Equation 1 using Lagrange multipliers, it is also unclear how large (or small) the value of alpha is and how it impacts the interpretations of the main result. \n\nc. Considering that one of the paper's main claims is a \"flat minima\" phenomenon, the paper would benefit from stronger theoretical results (apart from simulation-based arguments) surrounding this claim. \n\nd. Possible typo in main theorem. In the main theorem's statement for the nuclear norm, the result relies on $\\max(\\Sigma,\\alpha)$: should $\\Sigma$ be replaced with something like the maximum eigenvalue?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Na49L87kmP", "reviewer": "Reviewer_v15j", "summary": "This paper considers the trade-off between the variance (MSE) and bias (bias norm) for the regularized linear regression model, where the authors consider various Schatten norms for the regularized term. Then, the authors derive the expected MSE for the Gaussian and diagonal ensembles and compare it with the ridge regression model's one and show that the ridge regression model is not always the best in term of the trade-off by using some experiments on sythetic datasets.  There exit some similar works such as Bayatti and Montanari (2011), Samet et al. (2013) for Lasso, however the authors limit their work to the class of Schatten norms.", "questions": "How do your results in comparison with other norms which don't belong to Schatten class of norms such as Lasso (norm-$1$)?\n\nBesides, please check and correct the following typos and unprecise.\n\n+ $L \\in \\mathbb{R}^{k \\times N} \\rightarrow L \\in \\mathbb{R}^{d \\times N}$\n+ p.2, line 21 from the top: $\\mbox{var}\\_{\\epsilon}=L^T L$ should be changed to $\\mbox{var}\\_{\\epsilon}=\\sigma^2 L L^T$. To keep the later, you may change the definition of $\\hat{\\beta}$ to $L^T(X)Y$ throughout your paper.\n+ In the definition 1, you aim to minimize the variance subject to a constraint on bias by $C$. But, in the later (Theorem 2, Figure 1, etc.), it seems to me that you don't mention $C$ again, but only mention $\\alpha$. At least you should mention what is $\\alpha$ as a function of $C$ or vice versa in Theorem 2.\n+ In Figure 1, you only plot for the very special case when $X$ is diagonal. You should also plot for other cases of $X$. \n+ Typo in Theorem 2.1, the expectation should be over $\\epsilon$ only since you already take expectation over $X$ in MSE, and $Y$ is a function of $X$ and $\\epsilon$. \n+ Right below Figure 2: Figure 3.1 $\\rightarrow$ Figure 3. \n+ Similarly, Figure 3.2.1 $\\rightarrow$ Figure 5.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "+  The authors can obtain an exact expression for the average test error (MSE) for the spherical Gaussian model, and this bound is nearly matched to the experiment results (cf. Figure 2).\n+ Experiments show that Ridge regression, which used Frobenius norm, is not always the best option for the linear model in term of the trade-off between the variance (MSE) and the bias (the bias norm). More specifically, the authors show that the Frobenius norm and the nuclear one are likely to have the same average MSE on the spherical Gaussian ensemble (model) or the diagonal matrix model, but using the nuclear norm usually achieves better MSE than the Frobenius norm. This fact also holds when mapping the dataset via a random fourier transform (RFF).", "weaknesses": "+ The theoretical results only hold in the thermodynamic limit $N\\rightarrow \\infty, d \\rightarrow \\infty$ and $d/N\\rightarrow \\lambda$. This means that the results only hold when the number of observations is linear to the signal dimension. However, in common high-dimensional settings, the number of observations is usually sub-linear to the signal dimension. \n+ The result looks not an extension of the Gauss-Markov theorem since it only holds under expectation over $X$ when $X$ is an Gaussian ensemble or a diagonal ensemble. The Gauss-Markov theorem works for any $X$. \n+ It looks more interesting to compare your experiment results with other norms (outside the class of  Schatten norms) such as between the nuclear norm and Lasso (norm-$1$). \n+ Too many typos. Please check and correct them.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (poster)