# The Geometry of Over-parameterized Regression and Adversarial Perturbations

Jason W. Rocks
Department of Physics
Boston University
Boston, MA 02215
<EMAIL_ADDRESS>&Pankaj Mehta
Department of Physics
Faculty of Computing and Data Sciences
Boston University
Boston, MA 02215
<EMAIL_ADDRESS>Corresponding author.

###### Abstract

Classical regression has a simple geometric description in terms of a
projection of the training labels onto the column space of the design matrix.
However, for over-parameterized models – where the number of fit parameters is
large enough to perfectly fit the training data – this picture becomes
uninformative. Here, we present an alternative geometric interpretation of
regression that applies to both under- and over-parameterized models. Unlike
the classical picture which takes place in the space of training labels, our
new picture resides in the space of input features. This new feature-based
perspective provides a natural geometric interpretation of the double-descent
phenomenon in the context of bias and variance, explaining why it can occur
even in the absence of label noise. Furthermore, we show that adversarial
perturbations – small perturbations to the input features that result in large
changes in label values – are a generic feature of biased models, arising from
the underlying geometry. We demonstrate these ideas by analyzing three minimal
models for over-parameterized linear least squares regression: without basis
functions (input features equal model features) and with linear or nonlinear
basis functions (two-layer neural networks with linear or nonlinear activation
functions, respectively).

##  1 Introduction

Classical statistics has a long-standing practice of using geometry as a tool
to understand the behavior of linear regression models. Ordinary least squares
(OLS) regression has a particularly elegant interpretation: given a training
data set consisting of input features $\vec{\mathbf{x}}$ paired with labels
(responses) $y$, the fitted values $\hat{y}$ for the labels can be calculated
by orthogonally projecting them onto the column space of the design matrix [1]
(see Sec. 3 and Fig. 1). Geometric descriptions have also proven useful for
understanding related methods such as instrumental variables analysis, where
regression has a natural interpretation in terms of oblique projection [2].

In the classical setting, a model’s performance on new data points (the
generalization or test error) is well-understood in terms of the bias-variance
trade-off: while increasing a model’s complexity (e.g., increasing the number
of fit parameters) reduces bias (error due to the inability to fully express
patterns hidden in the data), it increases variance (error arising from over-
sensitivity to non-general aspects of the training set like noise). This
trade-off is reflected in the test error in the form of a classical,
“U-shaped” curve in which the test error first decreases with model complexity
until it reaches a minimum before increasing dramatically as the model
overfits the training data. Within this framework, optimal performance is
achieved at intermediate model complexities which strike a balance between
bias and variance [3, 4].

However, modern supervised learning techniques like Deep Learning methods defy
this dogma, achieving state-of-the-art performance using “over-parameterized”
models, where the number of fit parameters is large enough to perfectly fit
the training data [5]. In fact, increasing model complexity beyond the
interpolation threshold – the point at which the training error reaches zero –
actually results in a decrease in test error due to a drop in both bias and
variance, giving rise to what is now commonly referred to as “double-descent”
curves [6, 7] [see Fig. 1(a)]. Recently, the double-descent phenomenon has
been shown to be a fundamental property of over-parameterized models ranging
from modern neural networks to classical linear regression [8, 9, 10, 11, 12,
13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31,
32, 33, 34, 35, 36].

In this modern setting, the classical geometric picture suffers from two
significant drawbacks. First, the description is uninformative for over-
parameterized models since by definition, they always achieve zero training
error and hence the relevant projector simply becomes the identity. Second,
the classical picture fails to provide any intuition about test error and the
double-descent phenomenon.

In light of these inadequacies, here we present an alternative geometric
picture of OLS that resides in the space of input features rather than the
space of training labels. This picture applies to both under- and over-
parameterized models and provides simple geometric explanations for the
double-descent phenomenon and adversarial perturbations. We illustrate these
ideas by considering three different choices of basis functions for OLS
regression: (i) no basis functions (model features equal input features), (ii)
random linear features (a two-layer neural network with a linear activation
function and a random, fixed middle layer), and (iii) random nonlinear
features (a two-layer neural network with an arbitrary nonlinear activation
function and a random, fixed middle layer).

#### Summary of Results

We briefly summarize the contributions of this work:

  * •

We show that when an OLS model creates a label prediction $\hat{y}$ for a data
point $(y,\vec{\mathbf{x}})$, it also produces an internal representation of
the input features $\hat{\mathbf{x}}$. As a result, every data point is
associated with a pair of estimates for both the dependent and independent
variables $(\hat{y},\hat{\mathbf{x}})$.

  * •

We show that the internal representation of the input features can be
expressed as ${\hat{\mathbf{x}}=P_{f}\vec{\mathbf{x}}}$, where the matrix
operator $P_{f}$ first maps each data point onto the subspace
$\mathcal{F}_{W}$ of input features the model can express, and then onto the
subspace $\mathcal{F}_{X}$ of input feature spanned by the training data.

  * •

We discuss how the choice of basis functions for the model features affects
the geometric nature of the operator $P_{f}$ in the space of input features:
(i) a lack of basis functions results in an orthogonal projection, (ii) a
random linear basis results in an oblique projection, and (iii) a random
nonlinear basis results in what we call a “noisy” oblique projection.

  * •

We show how model prediction accuracy can be characterized by the angles
between directions in the subspaces $\mathcal{F}_{W}$ and $\mathcal{F}_{X}$
related by $P_{f}$. In particular, when one of these angles approaches
$90^{\circ}$, a model loses its ability to express a direction in the training
data, resulting in overfitting. We show how this behavior gives rise to the
double-descent phenomenon even in the absence of label noise.

  * •

We study the geometry of adversarial perturbations, showing that perturbations
to the input features of a data point decompose along “adversarial directions”
(the complement of the kernel $P_{f}$) that result in large typical changes in
the label prediction and “invariant directions” (the kernel $P_{f}$) that
result in no change on average. We show that the subspace of invariant
directions stems from model bias.

##  2 Theoretical Setup

In this work, we focus on the supervised learning task of OLS regression, in
which a linear model uses the relationships learned from a training data set
to accurately predict the labels of new data points based on their input
features. We consider a training data set of $M$ data points,
$\mathcal{D}=\\{(y_{a},\vec{\mathbf{x}}_{a})\\}_{a=1}^{M}$, each consisting of
a label $y$ and a vector of $N_{f}$ input features $\vec{\mathbf{x}}$. Each
label is related to its corresponding input features via the relationship
(teacher model)

$y(\vec{\mathbf{x}})=y^{*}(\vec{\mathbf{x}})+\varepsilon,$ (1)

where $y^{*}(\vec{\mathbf{x}})$ is an unknown function and $\varepsilon$ is
the label noise. For convenience, we organize the vectors of input features in
the training data set into the rows of a design matrix $X$ of size $M\times
N_{f}$ and define the length $M$ vectors of training labels $\vec{\mathbf{y}}$
and label noise $\vec{\boldsymbol{\varepsilon}}$.

Given an arbitrary data point $\vec{\mathbf{x}}$, the linear model for the
label predictions (student model) is then

$\hat{y}=\vec{\mathbf{z}}(\vec{\mathbf{x}})\cdot\hat{\mathbf{w}},$ (2)

where $\hat{\mathbf{w}}$ is a vector of $N_{p}$ fit parameters determined by
the training set. The vector of model features
$\vec{\mathbf{z}}(\vec{\mathbf{x}})$ is a deterministic function that
transforms the input features $\vec{\mathbf{x}}$ into a vector of $N_{p}$
basis functions (referred to as “hidden” features in the context of a two-
layer neural network). In analogy to the observation matrix, we organize the
vectors of model features evaluated on the input features of the training set,
$\\{\vec{\mathbf{z}}(\vec{\mathbf{x}}_{a})\\}_{a=1}^{M}$, into the rows of a
model feature matrix $Z$ of size $M\times N_{p}$ and define the length $M$
vector of training label predictions $\hat{\mathbf{y}}$.

We uniquely determine the fit parameters $\hat{\mathbf{w}}$ by requiring that
they minimize the following loss function evaluated on the training set
$\mathcal{D}$:

$L(\hat{\mathbf{w}};\mathcal{D})=\frac{1}{2}\norm{\Delta\vec{\mathbf{y}}}^{2}+\frac{\lambda}{2}\norm{\hat{\mathbf{w}}}^{2},$
(3)

where $\norm{\cdot}$ is the $L_{2}$ norm and
$\Delta\vec{\mathbf{y}}=\vec{\mathbf{y}}-\hat{\mathbf{y}}$ is the vector of
residual training label errors. The first term in the loss is the mean-squared
label error and the second term imposes $L_{2}$ regularization with
regularization parameter $\lambda$. To simplify the conceptual picture
presented in this work, we consider the ridge-less limit $\lambda\rightarrow
0$, allowing us to express the fit parameters that minimize Eq. (3) and the
resulting label predictions in Eq. (2), respectively, as

$\hat{\mathbf{w}}=Z^{+}\vec{\mathbf{y}},\qquad\hat{y}=\vec{\mathbf{z}}(\vec{\mathbf{x}})\cdot
Z^{+}\vec{\mathbf{y}},$ (4)

where + denotes the Moore-Penrose inverse, or pseudoinverse. Later, we discuss
the geometric interpretations of this solution for different sets of basis
functions for the model features $\vec{\mathbf{z}}(\vec{\mathbf{x}})$.

##  3 Geometry of Label Space

Figure 1: Double-descent phenomenon and corresponding geometric picture in
training label space $\mathcal{L}$. (a) The training error (blue squares) and
test error (black circles) for ordinary least square regression as a function
of the ratio of fit parameters $N_{p}$ to training data points $M$. In the
under-parameterized regime ($N_{p}<M$), increasing $N_{p}$ relative to $M$
decreases the training error until it reaches zero at the interpolation
threshold (vertical dashed line). At the same time, the test error follows a
classical “U-shaped” curve. In the over-parameterized regime ($N_{p}>M$), the
test error decreases while the training error remains zero. The $y$-axis is
scaled by the variance of the training labels $\sigma_{y}^{2}$. (b) In the
under-parameterized regime, the label predictions for the training set
$\hat{\mathbf{y}}$ are found by orthogonally projecting the training labels
$\vec{\mathbf{y}}$ onto the subspace $\mathcal{L}_{Z}\subset\mathcal{L}$
spanned by the rows of the matrix of model features evaluated on the training
set $Z$. (c) In the over-parameterized regime, the geometric picture is
trivial since the two spaces are identical, $\mathcal{L}_{Z}=\mathcal{L}$,
with the orthogonal projection resulting in label predictions with zero error,
$\hat{\mathbf{y}}=\vec{\mathbf{y}}$. See Sec. 5 and SI for numerical details.

Before introducing the geometric interpretation of the OLS solution in the
input feature space, we first review the classical picture of this solution in
the training label space. First, we define the training label space
$\mathcal{L}$ as the $M$-dimensional vector space spanned by all possible
values of the vector of training labels $\vec{\mathbf{y}}$ (e.g., here we
consider $\mathcal{L}=\mathbb{R}^{M}$ for continuous labels distributed over
$\mathbb{R}$). Based on Eq. (4), the vector of predicted labels takes the form

$\hat{\mathbf{y}}=P_{\ell}\vec{\mathbf{y}},\qquad P_{\ell}=ZZ^{+},$ (5)

where the $M\times M$ matrix $P_{\ell}$ represents an orthogonal projection
(${P_{\ell}^{2}=P_{\ell}=P_{\ell}^{T}}$) onto the subspace spanned by the
training data via the model features, or the rows of $Z$ which we denote as
$\mathcal{L}_{Z}$. In Fig. 1(b), we depict an example of this projection in
the classical under-parameterized setting for a synthetic data set where we
have represented $\mathcal{L}_{Z}$ as a line (see Sec. 5 and SI for details).

Most importantly, this picture provides a natural geometric form for the
training error,

$\displaystyle\mathcal{E}_{\mathrm{train}}$
$\displaystyle=\frac{1}{M}\norm{\Delta\vec{\mathbf{y}}}^{2},\qquad\Delta\vec{\mathbf{y}}=(I_{M}-P_{\ell})\vec{\mathbf{y}},$
(6)

where $I_{M}$ is the $M\times M$ identity matrix. From this, we see that the
residual training label errors $\Delta\vec{\mathbf{y}}$ are nonzero when
$\mathcal{L}_{Z}$ spans only a fraction of the possible directions in the
training label space $\mathcal{L}$.

While this description provides intuition for the under-parameterized regime,
its utility diminishes for over-parameterized models. As the complexity (e.g.,
the number of fit parameters) of a model is increased, the training error
decreases as the dimension of $\mathcal{L}_{Z}$ increases. Once the number of
parameters is large enough, the training error reaches zero at the
interpolation threshold as shown in Fig. 1(a). Past this threshold in the
over-parameterized, $\mathcal{L}_{Z}$ spans all of $\mathcal{L}$ and
$P_{\ell}$ trivially becomes the identity operator $I_{M}$, and as shown in
Fig. 1(c), this leads to zero residual label error,
$\vec{\mathbf{y}}=\hat{\mathbf{y}}$. Furthermore, this geometric picture
provides no intuition for the double-descent profile of the test error in Fig.
1(a) for neither the under- nor over-parameterized regimes.

##  4 Geometry of Feature Space

Figure 2:  Geometric picture in input feature space $\mathcal{F}$. (a) _No
Basis Functions_. To predict the label of a new data point, the input features
$\vec{\mathbf{x}}$ are orthogonally projected onto the subspace
$\mathcal{F}_{X}$ spanned by the input features of the training data (row
space of $X$) to obtain $\hat{\mathbf{x}}$, the datapoint as “seen” from the
perspective the model [see Eq. (9)]. (b)  _(Random) Linear Features_. The
input features $\vec{\mathbf{x}}$ undergo an oblique projection to obtain
$\hat{\mathbf{x}}$ which depends on the subspace $\mathcal{F}_{X}$, along with
$\mathcal{F}_{W}$, the subspace of input features that the model can express
(column space of $W$). To perform the oblique projection, the model first
orthogonally projects onto $\mathcal{F}_{W}$ to obtain $\hat{\mathbf{x}}_{W}$
and then continues orthogonally from $\mathcal{F}_{W}$ until it intersects
$\mathcal{F}_{X}$ to obtain $\hat{\mathbf{x}}$. (c) _(Random) Nonlinear
Features_. The input features $\vec{\mathbf{x}}$ undergo a “noisy” oblique
projection to obtain $\hat{\mathbf{x}}$, which again depends on the two
subspaces $\mathcal{F}_{X}$ and $\mathcal{F}_{W}$. Similar to an oblique
projection, $\hat{\mathbf{x}}$ is first orthogonally projected onto
$\mathcal{F}_{W}$. However, unlike an oblique projection, the mapping from
$\mathcal{F}_{W}$ to $\mathcal{F}_{X}$ has an extra component parallel to
$\mathcal{F}_{W}$, introducing additional error. (d) Mapping between subspaces
performed by each term in the singular value decomposition of $P_{f}$. The
subspace orientation angle $\theta_{i}$ measures the relative orientation
between the left and right singular vectors, $\vec{\mathbf{f}}_{W,i}$ and
$\vec{\mathbf{f}}_{X,i}$, while the projection deviation angle
$\delta\phi_{i}$ is zero when the mapping is orthogonal to
$\vec{\mathbf{f}}_{W,i}$. See Sec. 4 for details.

Next, we introduce the geometric interpretation of the OLS solution in the
input feature space $\mathcal{F}$. We define this space as the
$N_{f}$-dimensional vector space spanned by all possible values the of input
features $\vec{\mathbf{x}}$ (e.g., here we consider
$\mathcal{F}=\mathbb{R}^{N_{f}}$ for continuous input features distributed
over $\mathbb{R}$). In the SI, we formalize all following results as theorems
with corresponding proofs when appropriate.

#### Geometric Decomposition of Label Predictions

To develop our geometric description, we first decompose the labels and model
features into their linear and nonlinear components. For a given data point
$(y,\vec{\mathbf{x}})$, we define the label decomposition

$y(\vec{\mathbf{x}})=\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}+\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}})+\varepsilon,\qquad\vec{\boldsymbol{\beta}}\equiv\Sigma_{\vec{\mathbf{x}}}^{-1}\mathrm{Cov}_{\vec{\mathbf{x}}}\quantity[\vec{\mathbf{x}},y(\vec{\mathbf{x}})].$
(7)

The first term of $y(\vec{\mathbf{x}})$ captures the linear correlations
between the input features $\vec{\mathbf{x}}$ and the labels $y^{*}$,
characterized by the vector of $N_{f}$ parameters $\vec{\boldsymbol{\beta}}$.
The notation $\mathrm{Cov}_{\vec{\mathbf{x}}}[\cdot,\cdot]$ represents the
covariance evaluated with respect to the distribution of input features and we
define
${\Sigma_{\vec{\mathbf{x}},ij}\equiv\mathrm{Cov}_{\vec{\mathbf{x}}}[x_{i},x_{j}]}$
as the covariance matrix of the input features (assumed to be invertible). The
second term $\delta y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}})$ captures the
remaining nonlinear component of $y^{*}$ [defined as ${\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}})\equiv
y^{*}(\vec{\mathbf{x}})-\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}}$] and
is statistically independent of the linear term with respect to the input
features (see SI for proof).

Next, we define the model feature decomposition

$\vec{\mathbf{z}}(\vec{\mathbf{x}})=W^{T}\vec{\mathbf{x}}+\delta\vec{\mathbf{z}}_{\mathrm{NL}}(\vec{\mathbf{x}}),\qquad
W\equiv\Sigma_{\vec{\mathbf{x}}}^{-1}\mathrm{Cov}_{\vec{\mathbf{x}}}\quantity[\vec{\mathbf{x}},\vec{\mathbf{z}}^{T}(\vec{\mathbf{x}})].$
(8)

Analogous to the label decomposition, the first term of
$\vec{\mathbf{z}}(\vec{\mathbf{x}})$ captures the linear correlations between
the input features $\vec{\mathbf{x}}$ and the model features
$\vec{\mathbf{z}}$, characterized by the $N_{f}\times N_{p}$ matrix $W$. The
second term $\delta\vec{\mathbf{z}}_{\mathrm{NL}}(\vec{\mathbf{x}})$ captures
the remaining nonlinear component of $\vec{\mathbf{z}}$ [defined as
${\delta\vec{\mathbf{z}}_{\mathrm{NL}}(\vec{\mathbf{x}})\equiv\vec{\mathbf{z}}(\vec{\mathbf{x}})-W^{T}\vec{\mathbf{x}}}$],
and is also statistically independent of the linear term with respect to the
input features (see SI for proof).

Using the definitions above, we decompose the solution for the predicted
labels in Eq. (4) for an arbitrary data point $\vec{\mathbf{x}}$ as follows
(see SI for proof):

$\begin{gathered}\hat{y}(\vec{\mathbf{x}})=\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}+\delta\hat{y}(\vec{\mathbf{x}}),\qquad\hat{\mathbf{x}}\equiv
P_{f}\vec{\mathbf{x}},\qquad P_{f}\equiv(WZ^{+}X)^{T}\\\
\delta\hat{y}(\vec{\mathbf{x}})\equiv\delta\vec{\mathbf{z}}_{\mathrm{NL}}^{T}(\vec{\mathbf{x}})Z^{+}\vec{\mathbf{y}}+\vec{\mathbf{x}}^{T}WZ^{+}\quantity(\delta\vec{\mathbf{y}}^{*}_{\mathrm{NL}}+\vec{\boldsymbol{\varepsilon}}),\end{gathered}$
(9)

where we have defined the length $N_{f}$ vector $\hat{\mathbf{x}}$ in terms of
the $N_{f}\times N_{f}$ matrix operator $P_{f}$. The quantity
$\delta\hat{y}(\vec{\mathbf{x}})$ captures the elements of the labels or model
features that the model interprets as “noise” and goes to zero in the absence
of label noise and nonlinearities.

#### Internal Representation of Input Features

The vector $\hat{\mathbf{x}}$ and the operator $P_{f}$ play central roles in
our geometric description. Notice that in the expressions above,
$\hat{\mathbf{x}}$ is the counterpart to the label prediction $\hat{y}$, so
that for a given data point $(y,\vec{\mathbf{x}})$, the model produces a pair
of estimates for both the dependent and independent variables
$(\hat{y},\hat{\mathbf{x}})$. For this reason, one can interpret
$\hat{\mathbf{x}}$ as the model’s internal representation of the input
features $\vec{\mathbf{x}}$. In other words, $\hat{\mathbf{x}}$ is the data
point as “seen” from the perspective of the model via the operator $P_{f}$,
which encodes the patterns a model has learned from the training set $X$ based
on its model features as characterized by $W$.

#### Mapping Between Subspaces

In order to construct the estimate of the input features $\hat{\mathbf{x}}$
for a new data point $\vec{\mathbf{x}}$, the operator $P_{f}$ performs a
mapping from the subspace $\mathcal{F}_{W}$ of input features the model can
express (the row space of $W$) to the subspace $\mathcal{F}_{X}$ of input
features spanned by the training data set (the column space of $X$). To
analyze this mapping, we represent $P_{f}$ in terms of its singular value
decomposition (SVD),

$P_{f}=\sum_{i:\sigma_{i}>0}\sigma_{i}\vec{\mathbf{f}}_{X,i}\vec{\mathbf{f}}_{W,i}^{T},$
(10)

where the sum ranges over the non-zero singular values $\sigma_{i}$. The
collections of left and right singular vectors, $\vec{\mathbf{f}}_{X,i}$ and
$\vec{\mathbf{f}}_{W,i}$, each form orthonormal bases within $\mathcal{F}_{X}$
and $\mathcal{F}_{W}$, respectively. However, we note that it is possible for
these two bases to span only a portion of their respective subspaces if they
differ in dimension. As a result, the vectors $\vec{\mathbf{f}}_{X,i}$ span
only a subspace of $\mathcal{F}_{X}$ if the model completely misses some of
the input features spanned by the training data set. Similarly, the vectors
$\vec{\mathbf{f}}_{W,i}$ span only a subspace of $\mathcal{F}_{W}$ if some of
the input features the model can express are not actually used to learn
anything from the training set.

To understand the behavior of $P_{f}$, we consider the mapping independently
performed by each term in its SVD, depicted in Fig. 2(d). In the first step,
$\vec{\mathbf{x}}$ is mapped onto $\mathcal{F}_{W}$, resulting in the
intermediate representation
${\hat{\mathbf{x}}_{W}=\sum_{i}\hat{\mathbf{x}}_{W,i}=\sum_{i}\vec{\mathbf{f}}_{W,i}\vec{\mathbf{f}}_{W,i}^{T}\vec{\mathbf{x}}}$.
Next, the data point is mapped from $\mathcal{F}_{W}$ to $\mathcal{F}_{X}$
resulting in
${\hat{\mathbf{x}}=\sum_{i}\hat{\mathbf{x}}_{i}=\sum_{i}\sigma_{i}\vec{\mathbf{f}}_{X,i}\vec{\mathbf{f}}_{W,i}^{T}\hat{\mathbf{x}}_{W,i}=\sum_{i}\sigma_{i}\vec{\mathbf{f}}_{X,i}\vec{\mathbf{f}}_{W,i}^{T}\vec{\mathbf{x}}}$.
Furthermore, we define the subspace orientation angle $\theta_{i}$ between
each pair of singular vectors, $\vec{\mathbf{f}}_{X,i}$ and
$\vec{\mathbf{f}}_{W,i}$, shown in Fig. 2(d). Collectively, these angles
characterize the relative orientation and overlap between the subspaces
$\mathcal{F}_{X}$ and $\mathcal{F}_{W}$. Later, we find that in special cases,
data points are mapped orthogonally from $\mathcal{F}_{W}$ to
$\mathcal{F}_{X}$. To capture this behavior, we define the projection
deviation angle $\delta\phi_{i}$, which measures how this mapping deviates
from $90^{\circ}$ for each term in the SVD. As shown in Fig. 2(d), we define
$\delta\phi_{i}$ as $90^{\circ}$ minus the angle between the vectors
$\vec{\mathbf{f}}_{W,i}$ and $\hat{\mathbf{x}}_{i}-\hat{\mathbf{x}}_{W,i}$
(for $\vec{\mathbf{x}}=\vec{\mathbf{f}}_{W,i}$).

Next, we examine the properties of the geometric operator $P_{f}$ for specific
choices of basis functions for the model features
$\vec{\mathbf{z}}(\vec{\mathbf{x}})$.

#### No Basis Functions: Orthogonal Projection

The simplest set of model features in linear regression are the input features
themselves (where $W=I_{N_{f}}$, the $N_{f}\times N_{f}$ identity matrix),

$\vec{\mathbf{z}}(\vec{\mathbf{x}})=\vec{\mathbf{x}},\qquad
P_{f}=X^{+}X=\sum_{i}\vec{\mathbf{f}}_{X,i}\vec{\mathbf{f}}_{X,i}^{T}.$ (11)

In this case, the operator $P_{f}$ reduces to an orthogonal projector in input
feature space (${P_{f}^{2}=P_{f}=P_{f}^{T}}$, with all $\theta_{i}=0$ and
$\delta\phi_{i}$ undefined, see SI for proof).

Because the model can express any input features
($\mathcal{F}_{W}=\mathcal{F}$), the accuracy of the model’s internal
representation $\hat{\mathbf{x}}$ is only limited by the input features
sampled by the training data set via $X$. In Fig. 2(a), we show how the
orthogonal projection takes place for a synthetic data set. The model
constructs $\hat{\mathbf{x}}$ by orthogonally projecting $\vec{\mathbf{x}}$
onto the subspace $\mathcal{F}_{X}$ of input features spanned by the training
data, allowing the model to compare the new data point to data it encountered
previously during training.

#### (Random) Linear Features: Oblique Projection

Next, we consider a basis in which the model features are linearly related to
the input features via a matrix $W$ [defined consistent with Eq. (8)],

$\vec{\mathbf{z}}(\vec{\mathbf{x}})=W^{T}\vec{\mathbf{x}},\qquad
P_{f}=(W(XW)^{+}X)^{T}=\sum_{i}\frac{1}{\cos\theta_{i}}\vec{\mathbf{f}}_{X,i}\vec{\mathbf{f}}_{W,i}^{T}.$
(12)

For concreteness, we choose $W$ to be a random matrix so that the model is
equivalent to a two-layer neural network with a linear activation function
(with a random, fixed middle layer). For this basis, $P_{f}$ becomes an
oblique projector, satisfying the condition for a projector,
${P_{f}^{2}=P_{f}}$, but not for orthogonality, ${P_{f}=P_{f}^{T}}$. The
singular vectors form mutually orthogonal bases such that
${\vec{\mathbf{f}}_{W,i}\cdot\vec{\mathbf{f}}_{X,j}=\delta_{ij}\cos\theta_{i}}$
with $\sigma_{i}=1/\cos\theta_{i}$ and $\delta\phi_{i}=0$ (see SI for proof).

Unlike the model features in the previous section, here the accuracy of the
internal representation $\hat{\mathbf{x}}$ for a data point $\vec{\mathbf{x}}$
is limited by both the input features sampled by the training set via $X$ and
the input features that the model can express via $W$. The resulting form of
$P_{f}$ as an oblique projector reflects this behavior, depending on both the
subspaces $\mathcal{F}_{X}$ and $\mathcal{F}_{W}$. In Fig. 2(b), we
demonstrate how the oblique projection takes place in two steps for a
synthetic data set (see Sec. 5 and SI). First, the input features
$\vec{\mathbf{x}}$ are orthogonally projected onto the subspace
$\mathcal{F}_{W}$, resulting in an intermediate representation
$\hat{\mathbf{x}}_{W}$ in terms of input features the model can express. Next,
$\hat{\mathbf{x}}_{W}$ is mapped onto the second subspace $\mathcal{F}_{X}$,
but in a manner orthogonal to the first subspace $\mathcal{F}_{W}$, allowing
the model to compare the current data point to those that it encountered
previously in training. As a result of this second step, information contained
in $\vec{\mathbf{x}}$ may be lost if the subspaces do not perfectly overlap.
This occurs if any pair of left and right singular vectors have a subspace
orientation angle $\theta_{i}>0$. In Sec. 5, we discuss how this misalignment
can cause large test errors at the interpolation threshold when one of the
angles $\theta_{i}$ approaches $90^{\circ}$ and the model loses the ability to
accurately express a direction in the training data.

#### (Random) Nonlinear Features: “Noisy” Oblique Projection

Lastly, we consider nonlinear basis functions in which the model features take
the form

$\vec{\mathbf{z}}(\vec{\mathbf{x}})=\varphi\quantity(W^{T}\vec{\mathbf{x}}),\qquad
P_{f}=(WZ^{+}X)^{T}=\sum_{i}\frac{\cos\delta\phi_{i}}{\cos(\theta_{i}+\delta\phi_{i})}\vec{\mathbf{f}}_{X,i}\vec{\mathbf{f}}_{W,i}^{T},$
(13)

where $\varphi$ is an arbitrary nonlinear activation function that is applied
to each matrix element and $W$ is a matrix. For example, ReLU activation
corresponds to $\varphi(a)=C\max(0,a)$ where $C$ is a normalization constant
defined so that Eq. (8) holds. If we choose $W$ to be random, the model is
equivalent to a two-layer neural network with a nonlinear activation function
(with a random, fixed middle layer). Unlike with the linear model features of
the previous two sections, $P_{f}$ is no longer a true projector as the
condition $P_{f}^{2}=P_{f}$ generally no longer holds. In addition, the left
and right singular vectors are no longer mutually orthogonal and the singular
values depend on both $\theta_{i}$ and $\delta\phi_{i}$ as shown above (see SI
for proof).

Here, the accuracy of the internal representation $\hat{\mathbf{x}}$ is again
limited by both the input features spanned by the training data set and the
input features that the model can express. In Fig. 2(c), we show how $P_{f}$
differs from the linear case for a synthetic data set with ReLU activation
(see Sec. 5 and SI). Similar to oblique projection, a new data point
$\vec{\mathbf{x}}$ is first orthogonally projected onto $\mathcal{F}_{W}$ to
obtain the intermediate representation $\hat{\mathbf{x}}_{W}$ in terms of
input features the model can express. However, the mapping from of
$\hat{\mathbf{x}}_{W}$ onto $\mathcal{F}_{X}$, used to compare the data point
$\vec{\mathbf{x}}$ to those encountered during training, is no longer
orthogonal to $\mathcal{F}_{W}$. As a result, even if the model is capable of
expressing all possible input features ($\mathcal{F}_{W}=\mathcal{F}$) and the
input feature space is fully sampled during training
($\mathcal{F}_{X}=\mathcal{F}$), the internal representation
$\hat{\mathbf{x}}$ may not attain perfect accuracy. Although $P_{f}$ is no
longer formally an oblique projector, we observe that it does share some
elements in common with one. In Sec. 5, we show that near the interpolation
threshold, $P_{f}$ approximates to an oblique projector, resulting in large
test errors at the interpolation threshold when one of the angles $\theta_{i}$
approaches $90^{\circ}$. Furthermore, it is known that nonlinear model
features in regression tend to behave like linear features with noise [30].
Thus, we refer to $P_{f}$ as a “noisy” oblique projector in this case.

##  5 Double-Descent in the Absence of Noise

Figure 3: Double-descent phenomenon and corresponding geometric picture in
input feature space. (a) The total test error (black circles), geometric test
error (green triangles), geometric bias (blue squares), and geometric variance
(red diamonds) for random nonlinear features as a function of the ratio of fit
parameters $N_{p}$ to training data points $M$. The $y$-axis is scaled by the
variance of the training labels $\sigma_{y}^{2}$. (b) The subspace orientation
angle $\theta_{\max}$ (black circles) and the projection deviation angle
$\delta\phi_{i}$ (blue squares) associated with the maximum singular value in
the SVD of $P_{f}$ for the model in (a). At the interpolation threshold
(vertical dashed line), $\theta_{\max}$ approaches $90^{\circ}$ while
$\delta\phi_{i}$ approaches $0^{\circ}$. (c) In the under-parameterized regime
($N_{p}<M$), $\theta_{\max}<90^{\circ}$, resulting in low geometric variance.
(d) Near the interpolation threshold ($N_{p}\approx M$), $\theta_{\max}\approx
90^{\circ}$, resulting in large geometric variance. (e) In the over-
parameterized regime ($N_{p}>M$), once again $\theta_{\max}<90^{\circ}$ with
low geometric variance. See Sec. 5 and SI for numerical details.

A surprising feature of double-descent in over-parameterized regression is
that the phenomenon can occur even in the absence of label noise and
nonlinearities in the labels or model features [i.e., $\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}})=0$ and $\varepsilon=0$ in Eq. (7) and
$\delta\vec{\mathbf{z}}_{\mathrm{NL}}(\vec{\mathbf{x}})=0$ in Eq. (8)] [34].
To show this, we use the decompositions of the labels in Eq. (7) and their
predictions in Eq. (9) to define the geometric test error as
$\mathcal{E}_{\mathrm{geom}}\equiv[\Delta\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]^{2}$,
which depends on the vector quantity
$\Delta\vec{\mathbf{x}}\equiv(I_{N_{f}}-P_{f})\vec{\mathbf{x}}$. This quantity
is analogous to the residual label error defined in Eq. (6), so we interpret
it as the residual error of the input feature estimate $\hat{\mathbf{x}}$.
Furthermore, we decompose the geometric test error into geometric bias and
geometric variance terms by averaging over the distribution of possible
training sets $\mathcal{D}$ (see SI for proof),

$\begin{gathered}\mathrm{E}_{\mathcal{D}}\quantity[\mathcal{E}_{\mathrm{geom}}]=\mathrm{Bias}^{2}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]+\mathrm{Var}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}],\\\
\mathrm{Bias}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]\equiv\mathrm{E}_{\mathcal{D}}[\Delta\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}],\qquad\mathrm{Var}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]\equiv\mathrm{E}_{\mathcal{D}}[(\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}})^{2}]-\mathrm{E}_{\mathcal{D}}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]^{2}.\end{gathered}$
(14)

In the absence of label noise and nonlinearities, the traditional definitions
[1] of the test error, bias, and variance reduce to these geometric versions
(see SI for proof).

In Fig. 3(a), we plot these three quantities, along with the total test error
${\mathcal{E}_{\mathrm{test}}=[y(\vec{\mathbf{x}})-\hat{y}(\vec{\mathbf{x}})]^{2}}$,
for a basis of random nonlinear features with ReLU activation
$\smash{\varphi(x)=2\max(0,x)}$ for a synthetic data set with linear labels
$\smash{y^{*}(\vec{\mathbf{x}})=\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}}$.
We sample the input features $\vec{\mathbf{x}}$ and label noise $\varepsilon$
for the training and test data sets (each of size $M=512$) i.i.d. from normal
distributions with mean zero and variances $\sigma_{X}^{2}/N_{f}$ and
$\sigma_{\varepsilon}^{2}$, respectively. For each simulation, we also sample
the parameters $\smash{\vec{\boldsymbol{\beta}}}$ and matrix $W$ i.i.d. from
normal distributions with mean zero and variances
$\smash{\sigma_{W}^{2}/N_{p}}$ and $\sigma_{\beta}^{2}$, respectively. We
choose a signal-to-noise ratio for the labels of
$\smash{\sigma_{X}^{2}\sigma_{\beta}^{2}/\sigma_{\varepsilon}^{2}=10}$ and a
regularization parameter of $\lambda=10^{-8}$. We use this general setup to
construct all numerical examples and schematics presented in this work. All
points in error curves are averaged over at least 500 independent simulations
with small error bars indicating the error on the mean. We also include more
comprehensive numerical results for all three models in the SI.

As expected, we find that the divergence in the geometric test error stems
from the geometric variance at the interpolation threshold. Furthermore, we
find that the over-parameterized regime performs better than its under-
parameterized counterpart due to a monotonic decrease in the geometric bias as
the number of fit parameters increases. From Eq. (14), we see that the bias
arises from residual input feature error. Just as $\Delta y$ decreases as
$P_{\ell}$ approaches the identity $I_{M}$, resulting in lower training error,
$\Delta\vec{\mathbf{x}}$ decreases when $P_{f}$ approaches the identity
$I_{N_{f}}$, resulting in lower bias (see SI for numerics).

The overfitting that gives rise to the double-descent behavior occurs when a
model is not able to accurately express directions in the training data using
its model features. To see this in our geometric picture, we make use of the
observation that the test error diverges at the interpolation threshold when
$Z$ contains a small singular value, or alternatively, $P_{f}$ contains a
large singular value [9]. Therefore, we focus on the maximum singular value
$\sigma_{\max}$ which diverges and dominates the SVD of $P_{f}$ in Eq. (10)
near the interpolation threshold (see SI for numerics). In Fig 3(b), we plot
the subspace orientation angle $\theta_{\max}$ and projection deviation angle
$\delta\phi_{\max}$ associated with $\sigma_{\max}$. We find that near the
interpolation threshold, $\theta_{\max}$ approaches $90^{\circ}$, while
$\delta\phi_{\max}$ approaches $0^{\circ}$, and as a result, $P_{f}$
approaches a true oblique projector (see SI for proof).

Figure 4: Adversarial perturbations. (a) A perturbation to the input features
of a data point can be decomposed into “adversarial” and “invariant”
directions, parallel and perpendicular, respectively, to the subspace of input
features that can be captured by the model. (b) The change in the predicted
and true labels for $200$ random adversarial (blue) and $200$ random invariant
(orange) perturbations applied to the same data point with correlations shown
as lines (see Sec. 5 and SI for numerical details).

Most importantly, we find that test error arises when directions in the
subspaces $\mathcal{F}_{W}$ and $\mathcal{F}_{X}$ related via $P_{f}$ become
misaligned $\theta_{\max}>0$. In Figs. 3(c)-(d), we illustrate this behavior
for the singular vectors and angles associated with $\sigma_{\max}$. When
corresponding directions in the subspaces $\mathcal{F}_{W}$ and
$\mathcal{F}_{X}$ do not perfectly overlap, the model is not able to
accurately express the training data. The closer $\theta_{\max}$ is to
$90^{\circ}$, the worse the model becomes as it tends to overfit the poorly
expressed directions in the training data, resulting in large geometric
variance and test errors.

##  6 Adversarial Perturbations

Our geometric picture also yields insights into the nature of adversarial
perturbations. Inspired by Ref. [37], we consider perturbing a data point,
$\vec{\mathbf{x}}\rightarrow\vec{\mathbf{x}}+\eta\hat{\mathbf{e}}$, where
$\hat{\mathbf{e}}$ is a random unit vector and $\eta$ measures the magnitude
of the perturbation. As depicted in Fig. 4(a), any perturbation can be
decomposed into components parallel and perpendicular to the kernel of
$P_{f}$. Mathematically, this is the statement that we can write
${\hat{\mathbf{e}}=\hat{\mathbf{e}}_{\parallel}+\hat{\mathbf{e}}_{\perp}}$
where
${\sum_{i}\vec{\mathbf{f}}_{W,i}\vec{\mathbf{f}}_{W,i}^{T}\hat{\mathbf{e}}_{\parallel}=\hat{\mathbf{e}}_{\parallel}}$
and
${\sum_{i}\vec{\mathbf{f}}_{W,i}\vec{\mathbf{f}}_{W,i}^{T}\hat{\mathbf{e}}_{\perp}=0}$.
Our picture suggests that perturbations
$\hat{\mathbf{e}}=\hat{\mathbf{e}}_{\parallel}$, should have a large effect on
the label predictions $\hat{y}$, whereas the model should be insensitive to
perturbations $\hat{\mathbf{e}}=\hat{\mathbf{e}}_{\perp}$, which result in
similar internal representations $\hat{\mathbf{x}}$ (with any variations
coming entirely from noise and nonlinear effects). For this reason, the
complement to the kernel of $P_{f}$ is the subspace of adversarial directions
and the kernel of $P_{f}$ is that of invariant directions. Furthermore, Eq.
(14) indicates that the latter subspace is a result of nonzero geometric bias.

To test this intuition, we calculate how the model predictions $\hat{y}$ and
true labels $y$ change in response to small random perturbation from both
subspaces. In Fig. 4(b), we plot the derivatives
$\derivative{\hat{y}(\vec{\mathbf{x}}+\eta\hat{\mathbf{e}})}{\eta}$ and
$\derivative{y(\vec{\mathbf{x}}+\eta\hat{\mathbf{e}})}{\eta}$ for random
perturbations that are purely adversarial or invariant (see SI for details).
We find that adversarial perturbations tend to significantly change the
predicted labels when the true labels change (correlation is shown as a blue
line), whereas invariant perturbations on average do not result in changes to
the predicted labels that correlate with the changes in true labels (orange
line).

##  7 Conclusions

In this paper, we presented an alternative geometric interpretation of
regression, residing in the space of input features so that it applies both
under- and over-parameterized models. This new picture provides a geometric
interpretation of the double-descent phenomenon and suggests a natural scheme
for decomposing perturbations into adversarial and invariant directions. It
will be interesting to see which, if any, of these geometric intuitions can be
extended to more complicated settings such as deep neural networks and
classification tasks.

## Broader Impacts

The authors do not believe this theoretical work will raise any ethical
concerns or will generate any adverse future societal consequences.

## Acknowledgments

This work was supported by NIH NIGMS grant 1R35GM119461 and a Simons
Investigator in the Mathematical Modeling of Living Systems (MMLS) award to
PM. The authors also acknowledge support from the Shared Computing Cluster
administered by Boston University Research Computing Services.

## References

  * [1] Christopher M. Bishop.  Pattern Recognition and Machine Learning.  Springer, 2006.
  * [2] Roger J. Bowden and Darrell A. Turkington.  Instrumental Variables.  Cambridge University Press, 1990.
  * [3] Yaser S. Abu-Mostafa, Malik Magdon-Ismail, and Hsuan-Tien Lin.  Learning From Data.  AMLBook, 2012.
  * [4] Trever Hastie, Robert Tibshirani, and Jermoe Friedman.  The Elements of Statistical Learning: Data Mining, Inference, and Prediction.  Springer, 2nd edition, 2016.
  * [5] Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, and Oriol Vinyals.  Understanding Deep Learning Requires Re-thinking Generalization.  International Conference on Learning Representations (ICLR), 2017\.
  * [6] Mikhail Belkin, Daniel Hsu, Siyuan Ma, and Soumik Mandal.  Reconciling modern machine-learning practice and the classical bias–variance trade-off.  Proceedings of the National Academy of Sciences, 116(32):15849–15854, 2019.
  * [7] Marco Loog, Tom Viering, Alexander Mey, Jesse H. Krijthe, and David M. J. Tax.  A brief prehistory of double descent.  Proceedings of the National Academy of Sciences, 117(20):10625–10626, 2020.
  * [8] Ben Adlam and Jeffrey Pennington.  Understanding double descent requires a fine-grained bias-variance decomposition.  Advances in Neural Information Processing Systems (NeurIPS), 33:11022–11032, 2020.
  * [9] Madhu S. Advani, Andrew M. Saxe, and Haim Sompolinsky.  High-dimensional dynamics of generalization error in neural networks.  Neural Networks, 132:428–446, 2020.
  * [10] Jimmy Ba, Murat Erdogdu, Taiji Suzuki, Denny Wu, and Tianzong Zhang.  Generalization of Two-layer Neural Networks: An Asymptotic Viewpoint.  International Conference on Learning Representations (ICLR), 2020\.
  * [11] Jean Barbier, Florent Krzakala, Nicolas Macris, Léo Miolane, and Lenka Zdeborová.  Optimal errors and phase transitions in high-dimensional generalized linear models.  Proceedings of the National Academy of Sciences, 116(12):5451–5460, 2019.
  * [12] Peter L. Bartlett, Philip M. Long, Gábor Lugosi, and Alexander Tsigler.  Benign overfitting in linear regression.  Proceedings of the National Academy of Sciences, 117(48):30063–30070, 2020.
  * [13] Mikhail Belkin, Daniel Hsu, and Ji Xu.  Two Models of Double Descent for Weak Features.  SIAM Journal on Mathematics of Data Science, 2(4):1167–1180, 2020\.
  * [14] Koby Bibas, Yaniv Fogel, and Meir Feder.  A New Look at an Old Problem: A Universal Learning Approach to Linear Regression.  IEEE International Symposium on Information Theory (ISIT), pages 2304–2308, 2019.
  * [15] Zeyu Deng, Abla Kammoun, and Christos Thrampoulidis.  A Model of Double Descent for High-Dimensional Logistic Regression.  IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 4267–4271, 2020.
  * [16] Stéphane D’Ascoli, Levent Sagun, and Giulio Biroli.  Triple descent and the two kinds of overfitting: Where & why do they appear?  arXiv:2006.03509, 2020.
  * [17] Stéphane D’Ascoli, Maria Refinetti, Giulio Biroli, and Florent Krzakala.  Double Trouble in Double Descent: Bias and Variance(s) in the Lazy Regime.  Proceedings of the 37th International Conference on Machine Learning (ICML), PMLR 119:2280–2290, 2020.
  * [18] Michał Dereziński, Feynman Liang, and Michael W. Mahoney.  Exact expressions for double descent and implicit regularization via surrogate random design.  Advances in Neural Information Processing Systems (NeurIPS), 33:5152–5164, 2020.
  * [19] Federica Gerace, Bruno Loureiro, Florent Krzakala, Marc Mézard, and Lenka Zdeborová.  Generalisation error in learning with random features and the hidden manifold model.  Proceedings of the 37th International Conference on Machine Learning (ICML), PMLR 119:3452–3462, 2020.
  * [20] Trevor Hastie, Andrea Montanari, Saharon Rosset, and Ryan J. Tibshirani.  Surprises in High-Dimensional Ridgeless Least Squares Interpolation.  arXiv:1903.08560, 2019.
  * [21] Arthur Jacot, Berfin Şimşek, Francesco Spadaro, Clément Hongler, and Franck Gabriel.  Implicit regularization of random feature models.  Proceedings of the 37th International Conference on Machine Learning (ICML), PMLR 119:4631–4640, 2020.
  * [22] Ganesh Ramachandra Kini and Christos Thrampoulidis.  Analytic Study of Double Descent in Binary Classification: The Impact of Loss.  IEEE International Symposium on Information Theory (ISIT), pages 2527–2532, 2020.
  * [23] Andrew K. Lampinen and Surya Ganguli.  An analytic theory of generalization dynamics and transfer learning in deep linear networks.  International Conference on Learning Representations (ICLR), 2019\.
  * [24] Zhu Li, Weijie J. Su, and Dino Sejdinovic.  Benign overfitting and noisy features.  arXiv:2008.02901, 2020.
  * [25] Tengyuan Liang and Alexander Rakhlin.  Just interpolate: Kernel “Ridgeless” regression can generalize.  Annals of Statistics, 48(3):1329–1347, 2020.
  * [26] Tengyuan Liang, Alexander Rakhlin, and Xiyu Zhai.  On the Multiple Descent of Minimum-Norm Interpolants and Restricted Lower Isometry of Kernels.  Proceedings of Thirty Third Conference on Learning Theory, PMLR 125:2683–2711, 2020.
  * [27] Zhenyu Liao, Romain Couillet, and Michael W Mahoney.  A random matrix analysis of random Fourier features: beyond the Gaussian kernel, a precise phase transition, and the corresponding double descent.  Advances in Neural Information Processing Systems (NeurIPS), 33, 2020.
  * [28] Licong Lin and Edgar Dobriban.  What causes the test error? Going beyond bias-variance via ANOVA.  arXiv:2010.05170, 2020.
  * [29] Partha P Mitra.  Understanding overfitting peaks in generalization error: Analytical risk curves for $l_{2}$ and $l_{1}$ penalized interpolation.  arXiv:1906.03667, 2019.
  * [30] Song Mei and Andrea Montanari.  The generalization error of random features regression: Precise asymptotics and double descent curve.  arXiv:1908.05355, 2019.
  * [31] Vidya Muthukumar, Kailas Vodrahalli, and Anant Sahai.  Harmless interpolation of noisy data in regression.  IEEE International Symposium on Information Theory (ISIT), pages 2299–2303, 2019.
  * [32] Preetum Nakkiran.  More Data Can Hurt for Linear Regression: Sample-wise Double Descent.  arXiv:1912.07242, 2019.
  * [33] Preetum Nakkiran, Gal Kaplun, Yamini Bansal, Tristan Yang, Boaz Barak, and Ilya Sutskever.  Deep Double Descent: Where Bigger Models and More Data Hurt.  International Conference on Learning Representations (ICLR), 2020\.
  * [34] Jason W. Rocks and Pankaj Mehta.  Memorizing without overfitting: Bias, variance, and interpolation in over-parameterized models.  arXiv:2010.13933, 2020.
  * [35] Ji Xu and Daniel Hsu.  On the number of variables to use in principal component regression.  Advances in Neural Information Processing Systems (NeurIPS), 32, 2019.
  * [36] Zitong Yang, Yaodong Yu, Chong You, Jacob Steinhardt, and Yi Ma.  Rethinking Bias-Variance Trade-off for Generalization of Neural Networks.  Proceedings of the 37th International Conference on Machine Learning (ICML), PMLR 119:10767–10777, 2020.
  * [37] Ian J. Goodfellow, Jonathon Shlens, and Christian Szegedy.  Explaining and harnessing adversarial examples.  International Conference on Learning Representations (ICLR), 2015\.

## Supporting Information:
The Geometry of Over-parameterized Regression and Adversarial Perturbations

##  S1 Mathematical Proofs

######  Theorem 1.

In the label decomposition of Eq. (7), the linear term
$\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}$ and the nonlinear term $\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}})$ are statistically independent with
respect to the distribution of input features $\vec{\mathbf{x}}$.

###### Proof.

We show the covariance of $\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}$ and
$\delta y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}})$ is zero:

$\displaystyle\mathrm{Cov}_{\vec{\mathbf{x}}}[\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}},\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}})]$
$\displaystyle=\mathrm{Cov}_{\vec{\mathbf{x}}}[\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}},y^{*}(\vec{\mathbf{x}})-\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]$
(S1.1)
$\displaystyle=\vec{\boldsymbol{\beta}}\cdot\mathrm{Cov}_{\vec{\mathbf{x}}}[\vec{\mathbf{x}},y^{*}(\vec{\mathbf{x}})]-\vec{\boldsymbol{\beta}}\cdot\mathrm{Cov}_{\vec{\mathbf{x}}}[\vec{\mathbf{x}},\vec{\mathbf{x}}^{T}]\vec{\boldsymbol{\beta}}$
$\displaystyle=\vec{\boldsymbol{\beta}}\cdot\Sigma_{\vec{\mathbf{x}}}\vec{\boldsymbol{\beta}}-\vec{\boldsymbol{\beta}}\cdot\Sigma_{\vec{\mathbf{x}}}\vec{\boldsymbol{\beta}}$
$\displaystyle=0.$

∎

######  Theorem 2.

In the model feature decomposition of Eq. (8), the linear term
$W^{T}\vec{\mathbf{x}}$ and the nonlinear term
$\delta\vec{\mathbf{z}}_{\mathrm{NL}}(\vec{\mathbf{x}})$ are statistically
independent with respect to the distribution of input features
$\vec{\mathbf{x}}$.

###### Proof.

We show the covariance of $W^{T}\vec{\mathbf{x}}$ and
$\delta\vec{\mathbf{z}}_{\mathrm{NL}}(\vec{\mathbf{x}})$ is zero:

$\displaystyle\mathrm{Cov}_{\vec{\mathbf{x}}}[W^{T}\vec{\mathbf{x}},\delta\vec{\mathbf{z}}_{\mathrm{NL}}(\vec{\mathbf{x}})^{T}]$
$\displaystyle=\mathrm{Cov}_{\vec{\mathbf{x}}}[W^{T}\vec{\mathbf{x}},\vec{\mathbf{z}}^{T}(\vec{\mathbf{x}})-\vec{\mathbf{x}}^{T}W]$
(S1.2)
$\displaystyle=W^{T}\mathrm{Cov}_{\vec{\mathbf{x}}}[\vec{\mathbf{x}},\vec{\mathbf{z}}^{T}(\vec{\mathbf{x}})]-W^{T}\mathrm{Cov}_{\vec{\mathbf{x}}}[\vec{\mathbf{x}},\vec{\mathbf{x}}^{T}]W$
$\displaystyle=W^{T}\Sigma_{\vec{\mathbf{x}}}W-W^{T}\Sigma_{\vec{\mathbf{x}}}W$
$\displaystyle=0.$

∎

######  Theorem 3.

The label predictions can be expressed in terms of the geometric decomposition
in Eq. (9).

###### Proof.

We start with the OLS solution for the label predictions in Eq. (4) and then
substitute the label decomposition in Eq. (7) and model feature decomposition
in Eq. (8) to get

$\displaystyle\hat{y}(\vec{\mathbf{x}})$
$\displaystyle=\vec{\mathbf{z}}(\vec{\mathbf{x}})\cdot Z^{+}\vec{\mathbf{y}}$
(S1.3)
$\displaystyle=(W^{T}\vec{\mathbf{x}}+\delta\vec{\mathbf{z}}_{\mathrm{NL}}(\vec{\mathbf{x}}))\cdot
Z^{+}(X\vec{\boldsymbol{\beta}}+\delta\vec{\mathbf{y}}^{*}_{\mathrm{NL}}+\vec{\boldsymbol{\varepsilon}})$
$\displaystyle=\vec{\mathbf{x}}\cdot
WZ^{+}X\vec{\boldsymbol{\beta}}+\delta\vec{\mathbf{z}}_{\mathrm{NL}}(\vec{\mathbf{x}})^{T}Z^{+}(X\vec{\boldsymbol{\beta}}+\delta\vec{\mathbf{y}}^{*}_{\mathrm{NL}}+\vec{\boldsymbol{\varepsilon}})+\vec{\mathbf{x}}^{T}WZ^{+}(\delta\vec{\mathbf{y}}^{*}_{\mathrm{NL}}+\vec{\boldsymbol{\varepsilon}})$
$\displaystyle=\vec{\mathbf{x}}\cdot
P_{f}^{T}\vec{\boldsymbol{\beta}}+\delta\vec{\mathbf{z}}_{\mathrm{NL}}(\vec{\mathbf{x}})^{T}Z^{+}\vec{\mathbf{y}}+\vec{\mathbf{x}}^{T}WZ^{+}(\delta\vec{\mathbf{y}}^{*}_{\mathrm{NL}}+\vec{\boldsymbol{\varepsilon}})$
$\displaystyle=\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}+\delta\hat{y}(\vec{\mathbf{x}}),$

where where have defined the vector $\hat{\mathbf{x}}\equiv
P_{f}\vec{\mathbf{x}}$, the matrix operator $P_{f}\equiv(WZ^{+}X)^{T}$, and
the scalar quantity
${\delta\hat{y}(\vec{\mathbf{x}})\equiv\delta\vec{\mathbf{z}}_{\mathrm{NL}}(\vec{\mathbf{x}})^{T}Z^{+}\vec{\mathbf{y}}+\vec{\mathbf{x}}^{T}WZ^{+}(\delta\vec{\mathbf{y}}^{*}_{\mathrm{NL}}+\vec{\boldsymbol{\varepsilon}})}$.
∎

######  Theorem 4.

The operator $P_{f}=X^{+}X$ for model features without basis functions is an
orthogonal projector whose singular value decomposition is

$P_{f}=\sum_{i}\vec{\mathbf{f}}_{X,i}\vec{\mathbf{f}}_{X,i}^{T}$ (S1.4)

with subspace orientation angles $\theta_{i}=0$ and undefined projection
deviation angles $\delta\phi_{i}$.

###### Proof.

We start with the general form of the singular value decomposition of $P_{f}$,
given in Eq. (10),

$P_{f}=\sum_{i}\sigma_{i}\vec{\mathbf{f}}_{X,i}\vec{\mathbf{f}}_{W,i}^{T},$
(S1.5)

where $\vec{\mathbf{f}}_{W,i}\cdot\vec{\mathbf{f}}_{W,j}=\delta_{ij}$,
$\vec{\mathbf{f}}_{X,i}\cdot\vec{\mathbf{f}}_{X,j}=\delta_{ij}$ and
$\sigma_{i}>0$.

It is evident that $P_{f}$ is an orthogonal projector since

$P_{f}^{2}=X^{+}XX^{+}X=X^{+}X=P_{f}$ (S1.6)

and

$P_{f}^{T}=(X^{+}X)^{T}=X^{+}X=P_{f}$ (S1.7)

using the properties of the pseudoinverse.

Combining these two properties,

$\displaystyle P_{f}$ $\displaystyle=P_{f}^{T}P_{f}$ (S1.8)
$\displaystyle\sum_{j}\sigma_{j}\vec{\mathbf{f}}_{X,j}\vec{\mathbf{f}}_{W,j}^{T}$
$\displaystyle=\sum_{jk}\sigma_{j}\sigma_{k}\vec{\mathbf{f}}_{W,j}\vec{\mathbf{f}}_{X,j}^{T}\vec{\mathbf{f}}_{X,k}\vec{\mathbf{f}}_{W,k}^{T}$
$\displaystyle\sum_{j}\sigma_{j}\vec{\mathbf{f}}_{X,j}\vec{\mathbf{f}}_{W,j}^{T}$
$\displaystyle=\sum_{k}\sigma_{k}^{2}\vec{\mathbf{f}}_{W,k}\vec{\mathbf{f}}_{W,k}^{T}.$

We then right-multiply this equation by the vector $\vec{\mathbf{f}}_{W,i}$ to
find

$\displaystyle\sigma_{i}\vec{\mathbf{f}}_{X,i}$
$\displaystyle=\sigma_{i}^{2}\vec{\mathbf{f}}_{W,i}$ (S1.9)
$\displaystyle\vec{\mathbf{f}}_{X,i}$
$\displaystyle=\sigma_{i}\vec{\mathbf{f}}_{W,i}.$

Since these two vector are unit vectors and $\sigma_{i}$ is positive, it must
be that $\vec{\mathbf{f}}_{X,i}=\vec{\mathbf{f}}_{W,i}$ and $\sigma_{i}=1$,
resulting in the singular value decomposition

$P_{f}=\sum_{i}\vec{\mathbf{f}}_{X,i}\vec{\mathbf{f}}_{X,i}^{T}.$ (S1.10)

The subspace orientation angles $\theta_{i}$ are found using their definition,

$\cos\theta_{i}=\vec{\mathbf{f}}_{X,i}\cdot\vec{\mathbf{f}}_{W,i}=1,$ (S1.11)

indicating that $\theta_{i}=0$.

In addition, the projection deviation angles $\delta\phi_{i}$ are defined as

$\cos\quantity(\frac{\pi}{2}-\delta\phi_{i})=\frac{\vec{\mathbf{f}}_{W,i}\cdot(\hat{\mathbf{x}}_{i}-\hat{\mathbf{x}}_{W,i})}{\norm{\hat{\mathbf{x}}_{i}-\hat{\mathbf{x}}_{W,i}}}$
(S1.12)

with $\vec{\mathbf{x}}=\vec{\mathbf{f}}_{W,i}$. However, the vector in the
numerator has zero magnitude,

$\hat{\mathbf{x}}_{i}-\hat{\mathbf{x}}_{W,i}=\sigma_{i}\vec{\mathbf{f}}_{X,i}-\vec{\mathbf{f}}_{W,i}=0,$
(S1.13)

so the angles $\delta\phi_{i}$ are undefined. ∎

######  Theorem 5.

The operator $P_{f}=(W(XW)^{+}X)^{T}$ for linear features is an oblique
projector whose singular value decomposition is

$P_{f}=\sum_{i}\frac{1}{\cos\theta_{i}}\vec{\mathbf{f}}_{X,i}\vec{\mathbf{f}}_{W,i}^{T}$
(S1.14)

with
$\vec{\mathbf{f}}_{W,i}\cdot\vec{\mathbf{f}}_{X,i}=\delta_{ij}\cos\theta_{i}$
and $\delta\phi_{i}=0$.

###### Proof.

We start with the general form of the singular value decomposition of $P_{f}$,
given in Eq. (10),

$P_{f}=\sum_{i}\sigma_{i}\vec{\mathbf{f}}_{X,i}\vec{\mathbf{f}}_{W,i}^{T},$
(S1.15)

where $\vec{\mathbf{f}}_{W,i}\cdot\vec{\mathbf{f}}_{W,j}=\delta_{ij}$,
$\vec{\mathbf{f}}_{X,i}\cdot\vec{\mathbf{f}}_{X,j}=\delta_{ij}$ and
$\sigma_{i}>0$.

Clearly, $P_{f}$ is an oblique projector since

$\displaystyle P_{f}^{2}$ $\displaystyle=(W(XW)^{+}X)^{T}(W(XW)^{+}X)^{T}$
(S1.16) $\displaystyle=X^{T}(W^{T}X^{T})^{+}W^{T}X^{T}(W^{T}X^{T})^{+}W^{T}$
$\displaystyle=X^{T}(W^{T}X^{T})^{+}W^{T}$ $\displaystyle=(W(XW)^{+}X)^{T}$
$\displaystyle=P_{f}.$

Using this property, we find

$\displaystyle P_{f}$ $\displaystyle=P_{f}^{2}$ (S1.17)
$\displaystyle\sum_{k}\sigma_{k}\vec{\mathbf{f}}_{X,k}\vec{\mathbf{f}}_{W,k}^{T}$
$\displaystyle=\sum_{kl}\sigma_{k}\sigma_{l}\vec{\mathbf{f}}_{X,k}\vec{\mathbf{f}}_{W,k}^{T}\vec{\mathbf{f}}_{X,l}\vec{\mathbf{f}}_{W,l}^{T}.$

We then left-multiply by $\vec{\mathbf{f}}_{X,i}^{T}$ and right-multiply by
$\vec{\mathbf{f}}_{W,j}$, giving us

$\displaystyle\sigma_{i}\delta_{ij}$
$\displaystyle=\sigma_{i}^{2}\vec{\mathbf{f}}_{W,i}^{T}\vec{\mathbf{f}}_{X,j},$
(S1.18)

from which we find that
$\vec{\mathbf{f}}_{W,i}\cdot\vec{\mathbf{f}}_{X,j}=\frac{1}{\sigma_{i}}\delta_{ij}$,
where we have used the fact that $\sigma_{i}$ is positive to choose the sign.

Based on the definition of the subspace orientation angles $\theta_{i}$,

$\cos\theta_{i}=\vec{\mathbf{f}}_{X,i}\cdot\vec{\mathbf{f}}_{W,i}=\frac{1}{\sigma_{i}},$
(S1.19)

so the singular value decomposition is

$P_{f}=\sum_{i}\frac{1}{\cos\theta_{i}}\vec{\mathbf{f}}_{X,i}\vec{\mathbf{f}}_{W,i}^{T}$
(S1.20)

with
$\vec{\mathbf{f}}_{W,i}\cdot\vec{\mathbf{f}}_{X,j}=\cos\theta_{i}\delta_{ij}$.

We split the analysis of the projection deviation angles $\delta\phi_{i}$ into
two cases. First, when $\theta_{i}>0$, we have

$\displaystyle\cos\quantity(\frac{\pi}{2}-\delta\phi_{i})$
$\displaystyle=\frac{\vec{\mathbf{f}}_{W,i}\cdot(\hat{\mathbf{x}}_{i}-\hat{\mathbf{x}}_{W,i})}{\norm{\hat{\mathbf{x}}_{i}-\hat{\mathbf{x}}_{W,i}}}$
(S1.21)
$\displaystyle=\frac{\vec{\mathbf{f}}_{W,i}\cdot(\sigma_{i}\vec{\mathbf{f}}_{X,i}-\vec{\mathbf{f}}_{W,i})}{\norm*{\sigma_{i}\vec{\mathbf{f}}_{X,i}-\vec{\mathbf{f}}_{W,i}}}$
$\displaystyle=\frac{1-1}{\norm*{\sigma_{i}\vec{\mathbf{f}}_{X,i}-\vec{\mathbf{f}}_{W,i}}}$
$\displaystyle=0.$

∎

However, when $\theta_{i}=0$, the denominator
$\norm{\hat{\mathbf{x}}_{i}-\hat{\mathbf{x}}_{W,i}}$ is zero, so
$\delta\phi_{i}$ is undefined. In this case, we are free to take
$\delta\phi_{i}$ to be zero.

######  Theorem 6.

The operator $P_{f}=(WZ^{+}X)^{T}$ for arbitrary nonlinear features has a
singular value decomposition of the form

$P_{f}=\sum_{i}\frac{\cos\delta\phi}{\cos(\theta_{i}+\delta\phi_{i})}\vec{\mathbf{f}}_{X,i}\vec{\mathbf{f}}_{W,i}^{T}.$
(S1.22)

###### Proof.

We start with the definition of the projection deviation angles
$\delta\phi_{i}$,

$\displaystyle\cos\quantity(\frac{\pi}{2}-\delta\phi_{i})$
$\displaystyle=\frac{\vec{\mathbf{f}}_{W,i}\cdot(\hat{\mathbf{x}}_{i}-\hat{\mathbf{x}}_{W,i})}{\norm{\hat{\mathbf{x}}_{i}-\hat{\mathbf{x}}_{W,i}}}$
(S1.23)
$\displaystyle=\frac{\vec{\mathbf{f}}_{W,i}\cdot(\sigma_{i}\vec{\mathbf{f}}_{X,i}-\vec{\mathbf{f}}_{W,i})}{\norm*{\sigma_{i}\vec{\mathbf{f}}_{X,i}-\vec{\mathbf{f}}_{W,i}}}$
$\displaystyle=\frac{\sigma_{i}\cos\theta_{i}-1}{\sqrt{1+\sigma_{i}^{2}-2\sigma_{i}\cos\theta_{i}}}$

where we have also used the definition of the subspace orientation angle
$\theta_{i}$ as
$\cos\theta_{i}=\vec{\mathbf{f}}_{X,i}\cdot\vec{\mathbf{f}}_{W,i}$. Next, we
solve the above equation for $\sigma_{i}$ to find

$\displaystyle\sigma_{i}$
$\displaystyle=\frac{\sin\quantity(\frac{\pi}{2}-\delta\phi_{i})}{\sin\quantity(\theta_{i}+\frac{\pi}{2}-\delta\phi_{i})}$
(S1.24)
$\displaystyle=\frac{\cos\delta\phi_{i}}{\cos(\theta_{i}+\delta\phi_{i})},$

where we have chosen the sign so that singular values diverge when
$\theta_{i}+\delta\phi_{i}=\frac{\pi}{2}$.

The singular value decomposition is then

$P_{f}=\sum_{i}\frac{\cos\delta\phi}{\cos(\theta_{i}+\delta\phi_{i})}\vec{\mathbf{f}}_{X,i}\vec{\mathbf{f}}_{W,i}^{T}.$
(S1.25)

∎

######  Theorem 7.

The geometric test error can be decomposed into the geometric bias and
geometric variance as defined in Eq. (14).

###### Proof.

To derive the geometric bias-variance decomposition in Eq. (14), we start with
the definition of the geometric test error. We find

$\displaystyle\mathrm{E}_{\mathcal{D}}[\mathcal{E}_{\mathrm{geom}}]$
$\displaystyle=\mathrm{E}_{\mathcal{D}}[(\Delta\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}})^{2}]$
(S1.26)
$\displaystyle=\mathrm{E}_{\mathcal{D}}[(\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}-\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}})^{2}]$
$\displaystyle=\mathrm{E}_{\mathcal{D}}[(\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}-\mathrm{E}_{\mathcal{D}}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]+\mathrm{E}_{\mathcal{D}}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]-\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}})^{2}]$
$\displaystyle=(\mathrm{E}_{\mathcal{D}}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]-\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}})^{2}+\mathrm{E}_{\mathcal{D}}[(\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}-\mathrm{E}_{\mathcal{D}}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}])^{2}]$
$\displaystyle\qquad+2(\mathrm{E}_{\mathcal{D}}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]-\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}})\mathrm{E}_{\mathcal{D}}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}-\mathrm{E}_{\mathcal{D}}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]]$
$\displaystyle=(\mathrm{E}_{\mathcal{D}}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]-\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}})^{2}+\mathrm{E}_{\mathcal{D}}[(\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}})^{2}]-\mathrm{E}_{\mathcal{D}}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]^{2}$
$\displaystyle=\mathrm{Bias}^{2}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]+\mathrm{Var}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}],$

where the geometric bias and variance are defined as

$\displaystyle\mathrm{Bias}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]$
$\displaystyle=\mathrm{E}_{\mathcal{D}}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]-\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}},$
(S1.27)
$\displaystyle\mathrm{Var}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]$
$\displaystyle=\mathrm{E}_{\mathcal{D}}[(\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}})^{2}]-\mathrm{E}_{\mathcal{D}}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]^{2}.$

∎

######  Theorem 8.

The geometric test error, geometric bias and geometric variance reduce to the
standard definitions as presented in Ref. [1] in the absence of label noise
and nonlinearities.

###### Proof.

We substitute the decompositions for the true and predicted labels in Eqs. (7)
and (9) into the standard definitions of test error, bias and variance,
allowing us to write them in terms of their geometric counterparts. The total
test error averaged over the training data set $\mathcal{D}$ and the label
noise of the test data point $\varepsilon$ (assumed to have zero mean) is

$\displaystyle\mathrm{E}_{\mathcal{D},\varepsilon}[\mathcal{E}_{\mathrm{test}}]$
$\displaystyle=\mathrm{E}_{\mathcal{D},\varepsilon}[(y(\vec{\mathbf{x}})-\hat{y}(\vec{\mathbf{x}}))^{2}]$
(S1.28)
$\displaystyle=\mathrm{E}_{\mathcal{D},\varepsilon}[(\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}+\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}})+\varepsilon-\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}-\delta\hat{y}(\vec{\mathbf{x}}))^{2}]$
$\displaystyle=\mathrm{E}_{\mathcal{D},\varepsilon}[(\Delta\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}+\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}})+\varepsilon-\delta\hat{y}(\vec{\mathbf{x}}))^{2}]$
$\displaystyle=\mathrm{E}_{\mathcal{D}}[(\Delta\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}})^{2}]+\mathrm{E}_{\mathcal{D}}[(\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}})-\delta\hat{y}(\vec{\mathbf{x}}))^{2}]$
$\displaystyle\qquad+\mathrm{E}_{\varepsilon}[\varepsilon^{2}]+2\mathrm{E}_{\mathcal{D}}[(\Delta\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}})(\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}})-\delta\hat{y}(\vec{\mathbf{x}}))]$
$\displaystyle=\mathrm{E}_{\mathcal{D}}[\mathcal{E}_{\mathrm{geom}}]+\mathrm{E}_{\mathcal{D}}[(\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}})-\delta\hat{y}(\vec{\mathbf{x}}))^{2}]$
$\displaystyle\qquad+\mathrm{E}_{\varepsilon}[\varepsilon^{2}]+2\mathrm{E}_{\mathcal{D}}[(\Delta\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}})(\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}})-\delta\hat{y}(\vec{\mathbf{x}}))].$

The bias is

$\displaystyle\mathrm{Bias}[\hat{y}(\vec{\mathbf{x}})]$
$\displaystyle=\mathrm{E}_{\mathcal{D}}[\hat{y}(\vec{\mathbf{x}})]-y^{*}(\vec{\mathbf{x}})$
(S1.29)
$\displaystyle=\mathrm{E}_{\mathcal{D}}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}+\delta\hat{y}(\vec{\mathbf{x}})]-\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}-\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}})$
$\displaystyle=-\mathrm{E}_{\mathcal{D}}[\Delta\vec{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]+\mathrm{E}_{\mathcal{D}}[\delta\hat{y}(\vec{\mathbf{x}})]-\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}})$
$\displaystyle=-\mathrm{Bias}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]+\mathrm{E}_{\mathcal{D}}[\delta\hat{y}(\vec{\mathbf{x}})]-\delta
y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}}).$

The variance is

$\displaystyle\mathrm{Var}[\hat{y}(\vec{\mathbf{x}})]$
$\displaystyle=\mathrm{Var}_{\mathcal{D}}[\hat{y}(\vec{\mathbf{x}})]$ (S1.30)
$\displaystyle=\mathrm{Var}_{\mathcal{D}}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}+\delta\hat{y}(\vec{\mathbf{x}})]$
$\displaystyle=\mathrm{Var}_{\mathcal{D}}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]+\mathrm{Var}_{\mathcal{D}}[\delta\hat{y}(\vec{\mathbf{x}})]+2\mathrm{Cov}_{\mathcal{D}}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}},\delta\hat{y}(\vec{\mathbf{x}})].$

Clearly, these three quantities reduce to their geometric counterparts when
$\varepsilon$, $\delta y^{*}_{\mathrm{NL}}(\vec{\mathbf{x}})$ and
$\delta\vec{\mathbf{z}}_{\mathrm{NL}}(\vec{\mathbf{x}})$ are zero. ∎

######  Theorem 9.

When the maximum single singular value $\sigma_{\max}$ of $P_{f}$ diverges,
the subspace orientation angle $\theta_{\max}$ approaches $90^{\circ}$, and
the projection deviation angle $\delta\phi_{\max}$ approaches $0^{\circ}$,
$P_{f}$ is approximately an oblique projector for arbitrary nonlinear model
features.

###### Proof.

We start with the general form of the singular value decomposition of $P_{f}$,
given in Eq. (10),

$P_{f}=\sum_{i}\sigma_{i}\vec{\mathbf{f}}_{X,i}\vec{\mathbf{f}}_{W,i}^{T},$
(S1.31)

where $\vec{\mathbf{f}}_{W,i}\cdot\vec{\mathbf{f}}_{W,j}=\delta_{ij}$,
$\vec{\mathbf{f}}_{X,i}\cdot\vec{\mathbf{f}}_{X,j}=\delta_{ij}$ and
$\sigma_{i}>0$.

Now suppose a single singular value becomes much larger than the rest,
$\sigma_{\max}$. In this case, the singular value decomposition approximates
to

$P_{f}\approx\sigma_{\max}\vec{\mathbf{f}}_{X,\max}\vec{\mathbf{f}}_{W,\max}^{T}.$
(S1.32)

In addition if $\delta\phi_{\max}\approx 0^{\circ}$, then
$\sigma_{\max}\approx\frac{1}{\cos\theta_{\max}}$, so

$P_{f}\approx\frac{1}{\cos\theta_{\max}}\vec{\mathbf{f}}_{X,\max}\vec{\mathbf{f}}_{W,\max}^{T}.$
(S1.33)

Squaring this operator then results in

$\displaystyle P_{f}^{2}$
$\displaystyle=\frac{1}{\cos^{2}\theta_{\max}}\vec{\mathbf{f}}_{X,\max}\vec{\mathbf{f}}_{W,\max}^{T}\vec{\mathbf{f}}_{X,\max}\vec{\mathbf{f}}_{W,\max}^{T}$
(S1.34)
$\displaystyle=\frac{1}{\cos\theta_{\max}}\vec{\mathbf{f}}_{X,\max}\vec{\mathbf{f}}_{W,\max}^{T}$
$\displaystyle=P_{f},$

so $P_{f}$ is approximately an oblique projector. ∎

##  S2 Numerical Details

###  S2.1 Subspace Schematics

To generate the schematics of the training label space in Figs. 1(b)-(c) and
the input feature space in Figs. 2(a)-(c) and 4(a), we used the numerical
scheme detailed in Sec. 5. In each of these figures, we plot the relative
magnitudes of the $\vec{\mathbf{x}}$, $\hat{\mathbf{x}}_{W}$, and
$\hat{\mathbf{x}}$ and the angles between them to scale. We then used these
vectors to orient the lines representing the different subspaces. In each
figure, we only depict the portions of the subspaces that have a one-to-one
mapping via $P_{f}$, since these are the only directions that contribute to
$\hat{\mathbf{x}}$.

We note that in Figs. 2(b)-(c) and 4(a) the vector $\vec{\mathbf{x}}$,
$\hat{\mathbf{x}}_{W}$, and $\hat{\mathbf{x}}$ are drawn to lie in the same
plane for visual clarity, even though this is generally not the case. While
the angles between adjacent vectors are depicted to scale, the angles between
non-adjacent vectors are not drawn accurately.

In addition, the intermediate vector in Figs. 2(c) and 4(a), which lies
between the subspaces $\mathcal{F}_{X}$ and $\mathcal{F}_{W}$, generally does
not lie in the same place as $\hat{\mathbf{x}}_{W}$ and $\hat{\mathbf{x}}$. In
this case, the vector is an approximate 2$d$ representation found by finding
the component of $\hat{\mathbf{x}}-\hat{\mathbf{x}}_{W}$ that is perpendicular
to $\hat{\mathbf{x}}_{W}$ and lies in the same plane as both
$\hat{\mathbf{x}}_{W}$ and $\hat{\mathbf{x}}$.

In the schematics of the singular vectors in Figs. 2(c) and 3(c)-(d), the
relative magnitudes of the vectors $\hat{\mathbf{x}}_{W,i}$ and
$\hat{\mathbf{x}}_{i}$ and the angles $\theta_{i}$ and $\delta\phi_{i}$ are
drawn to scale. In each figure, we show the singular vectors and angles
corresponding to the largest singular value of $P_{f}$.

In Table S1, we report the model, the ratio of input features $N_{p}$ to
training data points, and the ratio of fit parameters $N_{p}$ to training data
points used to generate each schematic. In each case we use $M=512$ training
data points and round $N_{f}$ and $N_{p}$ down to the nearest integer.

Table S1: Model Details for Projector Schematics Figure | Model | Input Features $N_{f}/M$ | Fit Parameters $N_{p}/M$
---|---|---|---
1(b) | Random Nonlinear Features | $1/4$ | $1/4$
1(c) | Random Nonlinear Features | $1/4$ | $3$
2(a) | No Basis Functions | $6/5$ | $6/5$
2(b) | Random Linear Features | $6/5$ | $2$
2(c) | Random Nonlinear Features | $6/5$ | $2$
2(d) | Random Nonlinear Features | $1/4$ | $6/5$
3(c) | Random Nonlinear Features | $1/4$ | $5/6$
3(d) | Random Nonlinear Features | $1/4$ | $1$
3(e) | Random Nonlinear Features | $1/4$ | $6/5$
4(a) | Random Nonlinear Features | $6/5$ | $3$

###  S2.2 Bias-Variance Decompositions

Here, we describe our scheme to efficiently generate plots of the geometric
bias and geometric variance. First, we express the squared geometric bias and
geometric variance for a single test data point in terms of averages over two
independent training data sets as

$\displaystyle\mathrm{Bias}^{2}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]$
$\displaystyle=\mathrm{E}_{\mathcal{D}_{1},\mathcal{D}_{2}}[(\Delta\vec{\mathbf{x}}_{1}\cdot\vec{\boldsymbol{\beta}})(\Delta\vec{\mathbf{x}}_{2}\cdot\vec{\boldsymbol{\beta}})]$
(S2.1)
$\displaystyle\mathrm{Var}[\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}}]$
$\displaystyle=\mathrm{E}_{\mathcal{D}_{1}}[(\hat{\mathbf{x}}\cdot\vec{\boldsymbol{\beta}})^{2}]-\mathrm{E}_{\mathcal{D}_{1},\mathcal{D}_{2}}[(\vec{\mathbf{x}}_{1}\cdot\vec{\boldsymbol{\beta}})(\vec{\mathbf{x}}_{2}\cdot\vec{\boldsymbol{\beta}})]$

where $\hat{\mathbf{x}}_{1}$ and $\hat{\mathbf{x}}_{2}$ are the label
predictions from training on the different training data sets.

As suggested by these formulas, during each simulation, we independently
generate two training data sets of equal size, $\mathcal{D}_{1}$ and
$\mathcal{D}_{2}$, and perform regression separately on both, resulting in two
sets of fit parameters $\hat{\mathbf{w}}_{1}$ and $\hat{\mathbf{w}}_{2}$,
respectively. Using the results from the first training set, we calculate the
training error, test error, and geometric test error. We then evaluate the
geometric bias and geometric variance using the prescription above, but only
for the fixed pair of training data sets (and averaged over the test data
set). We repeat this process for many independent simulations and then average
all error quantities. The resulting geometric bias and geometric variance then
approximate to those in Eq. (S2.1), averaged over all random quantities.

###  S2.3 Adversarial Perturbations

To generate Fig. 4(b), we used the same model used to create the schematic in
Fig. 4(a). We then generated a random test data point $\vec{\mathbf{x}}$,
along with a collection of $200$ random pairs of adversarial and invariant
perturbations. To create each pair of perturbations, we randomly generated a
vector $\vec{\mathbf{e}}$ drawn form the same distribution as
$\vec{\mathbf{x}}$and decomposed it into its components parallel and
perpendicular to $\mathcal{F}_{W}$,

$\displaystyle\hat{\mathbf{e}}_{\parallel}$
$\displaystyle=\frac{\sum_{i}\vec{\mathbf{f}}_{W,i}\vec{\mathbf{f}}_{W,i}^{T}\vec{\mathbf{e}}}{\norm*{\sum_{i}\vec{\mathbf{f}}_{W,i}\vec{\mathbf{f}}_{W,i}^{T}\vec{\mathbf{e}}}}$
(S2.2) $\displaystyle\hat{\mathbf{e}}_{\perp}$
$\displaystyle=\frac{(I_{N_{f}}-\sum_{i}\vec{\mathbf{f}}_{W,i}\vec{\mathbf{f}}_{W,i}^{T})\vec{\mathbf{e}}}{\norm*{(I_{N_{f}}-\sum_{i}\vec{\mathbf{f}}_{W,i}\vec{\mathbf{f}}_{W,i}^{T})\vec{\mathbf{e}}}}.$

Next, we calculated the derivatives of the true and predicted labels for each
perturbation using finite difference,

$\displaystyle\derivative{y(\vec{\mathbf{x}}+\eta\hat{\mathbf{e}})}{\eta}$
$\displaystyle\approx\frac{y(\vec{\mathbf{x}}+\eta\hat{\mathbf{e}})-y(\vec{\mathbf{x}})}{\eta}$
(S2.3)
$\displaystyle\derivative{\hat{y}(\vec{\mathbf{x}}+\eta\hat{\mathbf{e}})}{\eta}$
$\displaystyle\approx\frac{\hat{y}(\vec{\mathbf{x}}+\eta\hat{\mathbf{e}})-\hat{y}(\vec{\mathbf{x}})}{\eta},$

with $\eta=10^{-2}$.

###  S2.4 Additional Numerical Results

Here, we provide more detailed numerical results for the three models we
consider in this work. We show results for the model without basis functions
in Fig. S5, with random linear features in Fig. S6, and random nonlinear
features in Fig. S7. For each model, we show the training error, test error,
geometric test error, geometric bias, and geometric variance. We also plots
the Frobenius ($L_{2,2}$) norms of the orthogonal complements of the data
projector $I_{M}-P_{\ell}$ and feature operator $I_{N_{f}}-P_{f}$, along with
the minimum singular value of $Z$, $\sigma_{Z,\min}$ and maximum singular
value of $P_{f}$, $\sigma_{\max}$. For the model with basis functions, we also
plot the subspace orientation angle $\theta_{\max}$ and projection deviation
angle $\delta\phi_{\max}$ corresponding to the maximum singular value of
$P_{f}$. Results are shown as a function of $N_{p}/M=N_{f}/M$ for the model
without basis functions and as a function of both $N_{p}/M$ and $N_{f}/M$ for
the two models with basis functions. For each model, black dashed lines
indicate phase transitions between different types of solutions [34]. Each
point or pixel is averaged over at least 100 independent simulations.

Figure S5: No Basis Functions. (a) Average training error (blue squares) and
test error (black circles). (b) Average geometric test error (black circles),
geometric bias (blue squares), and geometric variance (red diamonds). (c)
Average Frobenius matrix norm of the orthogonal complements of the data
projector $I_{M}-P_{\ell}$ (blue squares) and feature operator
$I_{N_{f}}-P_{f}$ (black circles). (d) Average minimum singular value of $Z$
(black circles) and maximum singular value of $P_{f}$ (blue squares). Vertical
dashed lines indicate the interpolation threshold.  Figure S6: Random Linear
Features. Average (a) training error, (b) test error, (e) geometric test
error, (f) geometric bias, and (g) geometric variance. Average Frobenius
matrix norm of the orthogonal complements of (c) the data projector
$I_{M}-P_{\ell}$ and (d) the feature operator $I_{N_{f}}-P_{f}$ . Average (h)
minimum singular value of $Z$ and (i) maximum singular value of $P_{f}$.
Average (j) subspace orientation angle $\theta_{\max}$ and (k) projection
deviation angle $\delta\phi_{\max}$ associated with the maximum singular value
of $P_{f}$. The angle $\delta\phi_{\max}$ is taken to be zero when
$\theta_{\max}=0$ and $\sigma_{\max}=1$. Vertical and horizontal lines
indicate the interpolation threshold. The diagonal dashed lines indicate the
boundary between the geometrically biased and unbiased regimes.  Figure S7:
Random Nonlinear Features. Average (a) training error, (b) test error, (e)
geometric test error, (f) geometric bias, and (g) geometric variance. Average
Frobenius matrix norm of the orthogonal complements of (c) the data projector
$I_{M}-P_{\ell}$ and (d) the feature operator $I_{N_{f}}-P_{f}$. Average (h)
minimum singular value of $Z$ and (i) maximum singular value of $P_{f}$.
Average (j) subspace orientation angle $\theta_{\max}$ and (k) projection
deviation angle $\delta\phi_{\max}$ associated with the maximum singular value
of $P_{f}$. Vertical lines indicate the interpolation threshold.