Spurious Correlations in High Dimensional Regression: The Roles of Regularization, Simplicity Bias and Over-Parameterization

We thank the reviewer for the remarkable care in reviewing our work and for the positive evaluation. We address concerns below.

Lack of clarity when citing (Han&Xu2023):

To address this and ease the comparison of our claims with the results in Han&Xu2023, we will add to the appendix a discussion about the notation in the proof of Theorem 4.3. Specifically, we will connect our definition of the Gaussian sequence model $\hat\theta^\rho$ with the definition in Equation (1.5) in Han&Xu2023, and our definition for the test function $\varphi$ (see our line 693) with their notation ($\textup{g}$) in their Theorem 2.3.

Other Comments:

• 102R, 146L, 136R: Thanks for pointing these typos out, we will fix them in the revision.

• 242L: While, in general, the conditional covariance depends on the particular value of $\bar x$, the definition of the Schur complement of the matrix $\Sigma$ does not rely on any specific instance of the random variable $x$. In the multivariate Gaussian case, it turns out that the conditional covariance is also independent of the particular instance $\bar x$, but we opted for leaving this notation at first as we did not consider it a trivial fact. If the reviewer finds this confusing, we are happy to remove this notation in the revision.

• The arguments of $f(\theta, x)$ are sometimes swapped: Thanks for spotting this, we will fix it.

• 625, 677, 683, 688: Thanks for pointing out these typos, we will fix them in the revision.

• 647: Thanks for noticing this. While it is true that the notation of $P_y$ also appears in the body in line 188R, in line 647 we are providing the proof for a statement prior to that part, and we opted for redundancy to avoid confusion.

• 690: Thanks for spotting this typo. Also, when we write that $\mathcal C(\cdot)$ is $\|| \Sigma \|| _{op}$-Lipschitz, we mean that the Lipschitz constant of $\mathcal C(\cdot)$ is upper bounded by the value of $\|| \Sigma \|| _{op}$. If the reviewer finds this notation confusing, we can elaborate more on this statement in the revision of the work.

• Hastie&al.2022: thanks for noticing the typo, we will fix it in the revision.

Questions for Authors:

• 146L: We thank the reviewer for pointing this out. While the implication holds for Gaussian distributions, in 146L we still did not state Assumption 4.1 and the wording at this point could be misleading. We will fix it in the revision.

• 662: We thank the reviewer for the question. Let us consider the inequality taken from the wikipedia page on Weyl’s inequality: $\lambda_{i+j-1}(A+B) \leq \lambda_{i}(A) + \lambda_{j}(B),$ where $A$ and $B$ are two $d \times d$ symmetric matrices and $\lambda_{j}(\cdot)$ denotes the $j$-th largest eigenvalue. Then, one could set $A = n \Sigma - Z^\top Z$ and $B = Z^\top Z$. Taking $i = 1$ and $j = d$ gives our Equation (32).

• 206R: We thank the reviewer for pointing this out. Indeed the sentence might lead to confusing interpretations. We propose to rephrase it as: “Thus, for large $d, n$, we can theoretically analyze $\mathcal C(\hat \theta_{\textup{LR}}(\lambda))$ via the deterministic quantity $\mathcal C^\Sigma(\lambda)$, which, as highlighted by Equation (12), depends on $\theta^*$, the covariance of the data $\Sigma$, and the regularization $\lambda$ via the parameter $\tau(\lambda)$ introduced in (13).”

We thank the reviewer for recognizing the rigor, clarity and importance of our results. We address concerns below.

Joint requirements $p = \omega(n)$ and $\log p = \Theta(\log n)$:

We note that the regime $p = \omega(n)$ and $\log p = \Theta(\log n)$ formally includes all scalings where $p \sim n^q$ for q > 1, as the latter is equivalent to $\log p = q \log n$. Thus, our assumptions include the faster scalings mentioned in the claims and evidence section of the review.

In case the reviewer referred to a polynomial regime also between $d$ and $n$ (if $n = \Omega(d^l)$, with $l \geq 2$, the RF model learns more than the linear component of the target as indicated in Hu&al.2024 – the paper referenced by the reviewer), it is indeed true that the higher order component of the features would behave qualitatively differently, making the RF model qualitatively different from a regularized linear regression. We have opted to focus on the proportional regime $n = \Theta(d)$ due to its popularity in the literature (see e.g. Mei&Montanari2020, Hastie&al.2020) and due to its closeness to standard datasets in deep learning.

Relation to work on linear regression with (possibly) noisy feature maps:

We thank the reviewer for bringing to our attention the paper by Loureiro et al. ("Learning curves of generic features maps for realistic datasets with a teacher-student model"), which has indeed a similar setting to ours. Using the notation $z = [x, y]$ as in our paper, they consider a teacher-student setting where the labels are defined as a function of the feature $x$ (see their Eq. (1.2)), while the estimator $\hat \theta$ is obtained via ERM using only the (correlated) features $y$ (see their Eq. (1.3)). Then, their work is focused on studying the training and generalization error of the model that has access only to the partial information. Our setting, instead, looks at the ERM on both features, where the model has direct access also to the core features $x$. Due to the similarity with their setting, we will mention this related work and remark the differences in the revision of the paper.

Setting where $z = x + y$ instead of $z = [x, y]$:

We thank the reviewer for the insightful comment. Let us consider the model

$z = x + y, g = x^\top \theta^* + \varepsilon.$

Then, we have

$\mathcal C = \text{Cov} \left( (\tilde x + y)^\top \hat \theta, x^\top \theta^* \right) = \hat \theta^\top \Sigma_{yx} \theta^*,$

and this quantity could be studied via the analysis in Han&Xu2023 as in our current setting, considering that the covariance of the data will take the form

$\Sigma_{zz} = \Sigma_{xx} + \Sigma_{yy} + \Sigma_{xy} + \Sigma_{yx}.$

In a nutshell, we expect the analysis for this setting to provide a qualitative behaviour similar to that unveiled in the current version of our work. In fact, the experiments on Color-MNIST (which does not strictly follow the model $z = [x, y]$, as the color overlaps with the core feature pattern as in the model $z = x + y$) suggest that our conclusions hold beyond the setting of orthogonal features. We remark that in the setting $z = [x, y]$ the optimal solution $\hat \theta = \theta^*$ gives $\mathcal C = 0$, while this is not necessarily the case in the setting $z = x + y$.

We will add a discussion on this point in the revision.

We thank the reviewer for the positive comments and the several interesting suggestions for extensions. We answer questions and address concerns below. We will incorporate the discussions in the revision.

Comparison with (Bombari et al., 2024):

Our work concerns the problem of spurious correlations, where a trained model is tested on a newly sampled data-point independent from training data, and we crucially use the independence between $\tilde x$, $x$, $y$ and $\hat \theta$ (see Proposition 4.2 and Theorem 4.3). In contrast, (Bombari et al., 2024) do not consider the problem of spurious correlations, but rather the setting where spurious features in the training set are memorized by an over-parameterized model. This is discussed in the first paragraph of their introduction, and it is quantitatively apparent in their definition of memorization in Equation (3.7), where the covariance is computed comparing the trained model evaluated on a spurious feature contained in the training set $y_i$ and the corresponding label $g_i$.

In other words, the setting of our work is related to robustness to distribution shift, while (Bombari et al., 2024) focus on a setting where the individual training data are memorized, raising potential privacy concerns. Thus, the two works look at qualitatively very different problems.

This difference is reflected in the proof strategies. While our work shares with (Bombari et al., 2024) an approach based on concentration of measure (and, consequently, also technical lemmas), the proof techniques are fundamentally different. Our work relies on the characterization of the ridge estimator $\hat \theta$ provided by Han&Xu2023 for linear regression, and it transfers the insights to random features via a point-wise equivalence principle. In contrast, the argument of (Bombari et al., 2024) is based on showing concentration of the auxiliary quantity $\mathcal F(z_i^s, z_i)$ that serves as a proxy to characterize the amount of memorization for an individual sample.

Setting might not capture the behaviour of deep neural networks:

While it is true that our analysis covers only high-dimensional regression, Figure 5 (left) shows a degree of similarity for shallow networks. For more complex and deep models, we point to the empirical results in Sagawa&al.2020, where higher penalty terms are shown to decrease test accuracy on ResNet50 and Bert. While it is hard to provide an exact predictive theory for deep models, we do believe our approach to capture important statistical aspects of the phenomenon of spurious correlations in more general settings than the one we precisely study.

Extension to feature learning:

Following Ba&al.2022 and Moniri&al.2023, one could extend our results to the setting where the target is not a linear function of the inputs and one step of gradient descent on the feature map improves the representation. We also note that our experiments with neural networks show concordance in the qualitative behavior of high-dimensional regression and 2-layer networks with both layers being trained.

Heavy tailed data:

The recent work by Adomaityte&al.2024 (“High-dimensional robust regression under heavy-tailed data: asymptotics and universality”) considers heavy tailed data in high-dimensional regression: the covariates are isotropic Gaussian with variance sampled from a distribution with heavy tails. Note that our problem setup requires a non-isotropic covariance, so one would have to first generalize their analysis accordingly. Then, a possible direction would be to investigate how different tail weights (between core and spurious features) favor learning of spurious correlations.

Different training objectives:

Empirically, the identified trade-off between regularization and spurious correlations has been verified in prior work (Sagawa&al.2020) looking at models trained on classification tasks.

Theoretically, work by Montanari&al.2023 (“The generalization error of max-margin linear classifiers”) and Deng&al.2020 (“A model of double descent for high-dimensional binary linear classification”) provides the asymptotics for the generalization error of max-margin linear classifiers, also in the setting where classification is performed on a set of random features. For classification, we could define $\mathcal C$ as in our Equation (3), taking the $\text{sign}$ of the output of the model (which for classification represents the prediction of the model at test time).

Equation (5.6) in Montanari&al.2023 provides a set of fixed point equations giving the limit deterministic value of the maximum margin and the prediction error (see their Theorem 3), also for non-isotropic covariates. Then, to extend our results to this setting, one approach could be to follow the strategy as in part c) of the proof technique in their Section 5.3, with the difference of computing $\mathcal C$ instead of the generalization error.

We thank the reviewer for the positive evaluation and helpful comments. We address concerns below.

Improving experimentation:

Following the reviewer’s suggestion, we will add the following experiments to the revision.

In https://ibb.co/21W8qbLJ, we consider Color-MNIST including all digits (rather than a subset). We train a 2-layer network on all classes with one-hot encoding and MSE loss. Odd (even) digits are red (blue) with probability $(1+\alpha)/2$, and blue (red) with probability $(1-\alpha)/2$. To compute $\mathcal C$, at test time we consider the parity of the logit with the highest value, with respect to the color of the image. The two figures correspond to two values of $\alpha$-s and follow a similar profile as Figure 5 (left), showing the same qualitative behaviour of $\mathcal L$ and $\mathcal C$ with respect to $\lambda$ for the full Color-MNIST dataset (i.e., in the multi-class setting).

In https://ibb.co/zWKBHk5k, we repeat the experiment in Figure 2 (right) for multiple values of $\alpha$, reporting $\mathcal L$ and $\mathcal C$ with respect to $\lambda$ for linear regression. The curves behave as expected: for any value of $\lambda$, as $\alpha$ decreases, $\mathcal C$ decreases. Furthermore, the (in-distribution) test loss decreases as $\alpha$ increases, in agreement with our discussion at lines 347-350 (left).

We note that the choice of considering a data split in predictive and spurious features is conceptually similar to Section 5.1 of “Invariant risk minimization” by Arjovsky et al. As for our CIFAR-10 implementation, our experiment is designed to verify our claims on the simplicity bias in a controllable setting. In fact, introducing a tunable amount of noise allows us to modify $\lambda_{\max}(\Sigma_{yy})$ and use it as an independent variable in Figure 4. Nevertheless, considering image backgrounds as spurious features was done in the seminal cited work by Xiao&al.2020. Besides, the theoretical results suggesting that higher values of the regularizer $\lambda$ can be associated with higher values of $\mathcal C$ are also supported by the numerical evidence in Table 1 of Sagawa2020a.

We thank the reviewer for mentioning the Corrupted CIFAR-10 dataset considered e.g. in “Avoiding spurious correlations via logit correction” by Liu et al. In https://ibb.co/TBg7Khwn, we train a 2-layer network and a random feature model on the classes ‘trucks’ and 'boats’, enforcing a correlation in the training set with respectively the textures ‘brightness’ and ‘glass_blur’, as in the available data for C-CIFAR-10 (here we use correlation $\alpha = 0.95$). In both figures, we see a mild increase in $\mathcal C$ as $\lambda$ initially increases, until the later decrease predicted by Proposition 5.1. The profiles are also qualitatively similar to the ones of Figure 5 (left).

More discussion in light of literature:

The main difference with (Qiu et al., 2024; Pezeshki et al., 2011) is that these works mainly focus on how spurious correlations evolve during training, while our work studies spurious correlations at convergence.

To further connect with (Qiu et al., 2024; Pezeshki et al., 2011), we briefly discuss the intuition coming from our results from a dynamical perspective. In Section 5, we argue that $\lambda_{\max}(\Sigma_{yy})$ is related to $\mathcal C$, suggesting a measure of the simplicity of the feature $y$. In linear regression, solving gradient flow ($d \theta = - \nabla_\theta \mathcal L(\theta)dt$) gives

$\theta(t) = (1 - e^{ - (X^\top X + n \lambda I ) t}) \hat \theta.$

Thus, the components of $\hat \theta$ aligned with the top eigenspaces of $X^\top X$ converge earlier than the others. Hence, if $X^\top X \sim n \Sigma$, it is natural to expect that spurious features are learned faster the easier they are, and that they would prevail with respect to the core features (according to our bound in Proposition 5.1).

Other comments/suggestions:

We thank the reviewer for the detailed comments, which will improve the clarity of the revision:

• We will explicitly note that the test loss is in distribution after Equation (2).
• We will clarify the usage of the notation $\lambda$, as well as of $y$ to denote spurious features and $g$ to denote labels.
• We will fix the typo in line 34L.
• We will add the suggested remark before Theorem 4.3.
• We will discuss the difference between the two different scaling regimes considered in Theorem 4.3 and Proposition 4.2.