# Asymptotic Characterisation of Robust Empirical Risk Minimisation
Performance
in the Presence of Outliers

Matteo Vilucchio  Information, Learning and Physics Laboratory, EPFL, 1015
Lausanne, Switzerland  Emanuele Troiani  Statistical Physics of Computation
Laboratory, EPFL, 1015 Lausanne, Switzerland  Vittorio Erba  Statistical
Physics of Computation Laboratory, EPFL, 1015 Lausanne, Switzerland  Florent
Krzakala  Information, Learning and Physics Laboratory, EPFL, 1015 Lausanne,
Switzerland

###### Abstract

We study robust linear regression in high-dimension, when both the dimension
$d$ and the number of data points $n$ diverge with a fixed ratio $\alpha=n/d$,
and study a data model that includes outliers. We provide exact asymptotics
for the performances of the empirical risk minimisation (ERM) using
$\ell_{2}$-regularised $\ell_{2}$, $\ell_{1}$, and Huber loss, which are the
standard approach to such problems. We focus on two metrics for the
performance: the generalisation error to similar datasets with outliers, and
the estimation error of the original, unpolluted function. Our results are
compared with the information theoretic Bayes-optimal estimation bound. For
the generalization error, we find that optimally-regularised ERM is
asymptotically consistent in the large sample complexity limit if one perform
a simple calibration, and compute the rates of convergence. For the estimation
error however, we show that due to a norm calibration mismatch, the
consistency of the estimator requires an oracle estimate of the optimal norm,
or the presence of a cross-validation set not corrupted by the outliers. We
examine in detail how performance depends on the loss function and on the
degree of outlier corruption in the training set and identify a region of
parameters where the optimal performance of the Huber loss is identical to
that of the $\ell_{2}$ loss, offering insights into the use cases of different
loss functions.

##  I Introduction

The problem of robust regression — i.e. linear regression with spurious
outliers in the training set — has been extensively studied in the classical
low-dimensional setting [1]. More recently, the high-dimensional setting has
also been considered mainly concerning the characterisation of the asymptotic
performance of robust estimators — see Section II for a review of related
works. In this paper, we provide new theoretical results on the problem in the
high-dimensional regime where the dimensionality of the problem $d$ is
comparable to the number of training samples $n$ available ($n,d\to\infty$
with fixed ratio $\alpha=n/d$), investigating in particular the phenomenology
arising due to the presence of outliers.

We build upon the recent progresses in high-dimensional statistics that
characterise the asymptotic performance of the Bayes Optimal (BO) and the
Empirical Risk Minimisation (ERM) estimators for general data models learned
with generalised linear models. In particular, we consider a probabilistic
model of training datasets — first introduced in [2] — that features a tunable
presence of outliers and we investigate the high-dimensional setting
performances of some standard regression estimators. We focus our attention on
empirical risk minimisers under the $\ell_{2}$-regularised $\ell_{2}$,
$\ell_{1}$ and Huber losses [3], and study how they asymptotically perform in
various outlier settings (e.g. small/large outlier percentage, small/large
outlier variance). We focus on two metrics for the performance, the
generalisation error to similar datasets with outliers, and the estimation
error of the original, noiseless function, and we compare to information
theoretic Bayes-optimal bounds.

For the probabilistic model of outliers we consider, eq. (1), we provide
closed-form expressions for the high-dimensional limit of the generalisation
error of ERM for the $\ell_{2}$, $\ell_{1}$ and Huber losses, and for the BO
estimator (for which we compute exact rates at large sample complexity
$\alpha$). Our main results are:

  * •

Optimally-regularised ERM is not always consistent. Based on the value of the
parameters characterising the outliers’ distribution, optimally-regularised
ERM estimation can fail to achieve BO performance for large sample complexity
$\alpha$, i.e. fail to achieve consistency. This happens for both the
generalisation error and for the estimation error, in different regions of the
parameter space. Notably, for the generalisation error there is a regime in
which optimally-regularised ERM with $\ell_{1}$ loss or the Huber loss (with
un-optimised scale parameter) is not consistent, while optimally-regularised
ERM with $\ell_{2}$ and optimally-scaled Huber loss is consistent. This
demonstrates the superiority of optimally-scaled Huber loss over the other
losses.

  * •

Calibration of ERM minimiser norm restores consistency. In the aforementioned
cases in which ERM is not consistent, the inconsistency is caused exclusively
by a mismatch between the norm of the ERM minimiser and that of the ground-
truth. We show that consistency for the generalization error, and optimal
rates, can be recovered by cross-validating on a test set to find the optimal
norm.

  * •

Non-consistency of the estimation error. Even when the cross-validated
estimator is consistent and achieves optimal rates for the generalization
error, it may remain inconsistent for the estimation error on the original,
unpolluted function, which is arguably a more interesting objective. While
this, too, can be calibrated, such a calibration requires either having access
to an oracle knowledge of the BO norm, or to a clean, unpolluted dataset to
cross-validate for the optimal norm. In the absence of such an ideal hold-out
dataset, the Huber and $\ell_{1}$ methods do not lead to consistent estimates
of the objective functions, and robust estimation is biased.

  * •

Dependency of the generalisation error on the parameters. We study the
dependency of the generalisation error on the amount and strength of the
outliers in the training dataset. We find that there is an extended region of
parameters where the performance of the Huber loss is the same as the simpler
$\ell_{2}$ function, offering insights into the use cases of different loss
functions.

In Section II we discuss related works. In Section III we introduce the data
model we consider, and the noise model accounting for outliers. In Section IV
we provide the analytical characterisation for the estimation and
generalisation errors for the ERM and BO estimators. In Section V we discuss
the large sample complexity limit, and prove our results concerning the
consistency of ERM estimators. In Section VI we discuss the dependency of the
generalisation error on the parameters characterising the outliers’
distribution. We provide more details and proofs figures in the Supplemental
Materials. The code to reproduce the figures can be found as a package
available at https://github.com/IdePHICS/RobustRegression.

##  II Related works

The problem of regression in the presence of outliers has been studied
extensively, starting with Huber’s seminal works [3, 4] where he first
introduced M-estimators as a robust alternative to classical estimators and
secondly extended the concept of robustness to hypothesis testing. A
comprehensive framework for understanding and defining robustness in
statistical procedures was then given by [5, 6].

Several algorithms have been proposed to approach problems with outliers — see
the review [1]. In most cases, these are variants of standard linear
regression based, for example, on custom loss functions [3], on hard-
thresholding [7], on some sparsity assumption [8] or on the sum of squares
method [9].

The choice of the outlier model is crucial for theoretical analysis and as an
assumption to propose practical algorithms. The double Gaussian noise model
that we consider has been first introduced in [2]. A similar model that
considers outliers as a two-way process is studied in [10]. Recently this
model has been reconsidered by [11], where they provide bounds on the
generalisation error in the case of Huber-loss regression. Other models have
also been proposed; see, for example, [12, 13, 14, 15, 16].

There have been numerous works analysing the generalisation error for
M-estimation and robust regression in high dimensions. In [17, 12, 18], they
consider high-dimensional M-estimation for a model of data with exponentially-
distributed outlier noise. They provide results on exact asymptotics and study
the trade-off between using $\ell_{2}$ and $\ell_{1}$ loss. In [19], they
prove bounds for the precise asymptotics of the generalisation error of
$\ell_{2}$ regularised $\ell_{2}$ loss, while in [20], they prove asymptotics
for M-estimation using Gordon min-max techniques. In [21], they consider a
sparse model under an almost-high-dimensional limit $n\sim\log(d)$. In [22], a
comparison of the precise asymptotics of $\ell_{1}$-regularised $\ell_{2}$,
$\ell_{1}$ and Huber loss is presented, with a focus on the interplay between
the sparsity and the sample complexity of the problem. In [23], they present
results on the asymptotic normality of M-estimators.

The characterisation of the generalisation error for ERM that we use and
specialise to our specific data model eq. (1) has been proven in [24] and more
recently in greater generality in [25]. The characterisation of the Bayes
optimal performance has been rigorously established in [26]. We are using
these results as the basis of our investigations in this paper. We stress that
all these works do not consider model outliers, thus providing no direct
insight on robust regression.

##  III Data and Noise Model

We will consider a linear regression task in $\mathbb{R}^{d}$, i.e. we will
fit a linear model — the “student” —
$\hat{y}(\boldsymbol{x};\boldsymbol{w})=\boldsymbol{x}\cdot\boldsymbol{w}/\sqrt{d}$
to a dataset $\mathcal{D}=\quantity{(\boldsymbol{x}_{i},y_{i})}_{i=1}^{n}$
made of $n$ pairs of data-points $\boldsymbol{x}_{i}\in\mathbb{R}^{d}$ and
labels $y_{i}\in\mathbb{R}$ sampled independently from a joint probability
distribution $P(\boldsymbol{x},y)$.

We will consider a probabilistic generative model $P(\boldsymbol{x},y)=P_{\rm
out}(y|\boldsymbol{x}\cdot\boldsymbol{w}^{\star}/\sqrt{d})P_{\rm
input}(\boldsymbol{x})$. In this model, the dataset labels are generated by
corrupting the output of a “teacher” linear model with weights
$\boldsymbol{w}^{*}$ through an output channel $P_{\rm out}$. We will consider
the simple Gaussian case where the teacher vector is distributed as
$\boldsymbol{w}^{\star}\sim\mathcal{N}(0,\mathbbm{1}_{d})$, and the samples
$\boldsymbol{x}$ have Gaussian distribution $P_{\rm
input}=\mathcal{N}(0,\mathbbm{1}_{d})$. In order to model datasets with
outliers, we consider the following generative process for the labels

$P_{\rm
out}(y|y^{\star})=\begin{cases}y^{\star}+z\sqrt{\Delta_{\text{\tiny{IN}}}}&\text{with
prob. }1-\epsilon\\\ \beta
y^{\star}+z\sqrt{\Delta_{\text{\tiny{OUT}}}}&\text{with prob. }\epsilon\\\
\end{cases}$ (1)

where $y^{\star}=\boldsymbol{x}\cdot\boldsymbol{w}^{\star}/\sqrt{d}$ is the
clean teacher-generated label, $z$ is a standard Gaussian noise,
$\Delta_{\text{\tiny{IN}}}$ is the variance of the non-outlying training
points, $\Delta_{\text{\tiny{OUT}}}$ is the variance of the outlying training
points, $\epsilon$ is the fraction of outliers and $\beta\geq 0$ is a
parameter controlling the outliers’ norm. In the following, we will consider
often the special case of $\beta=0$, but all of our results can be extended
for any value $\beta\in\mathbb{R}_{+}$ (see Appendix A.2).

We will consider two different metrics to analyse the performance of the
estimators. The first metric is the generalisation error of a model with
weights $\boldsymbol{w}$, defined as

$E_{\mathrm{gen}}(\boldsymbol{w})=\mathbb{E}_{(\boldsymbol{x}_{\text{new}},y_{\text{new}})}\quantity[\quantity(y_{\text{new}}-\frac{\boldsymbol{x}_{\text{new}}\cdot\boldsymbol{w}}{\sqrt{d}})^{2}]\,,$
(2)

where the average is over a sample-label pair
$(\boldsymbol{x}_{\text{new}},y_{\text{new}})$ generated with the same
distribution as the training set. The second metric we will consider is the
estimation error, which is the mean-square error between the teacher and
student weights, defined as

$E_{\mathrm{estim}}(\boldsymbol{w})=\frac{1}{d}\norm{\boldsymbol{w}^{\star}-\boldsymbol{w}}_{2}^{2}\,,$
(3)

We remark that, as the sample distribution $P_{\rm input}(\boldsymbol{x})$ is
Gaussian, the estimation error also coincides with the generalisation error
evaluated on a clean dataset

$E_{\mathrm{estim}}(\boldsymbol{w})=\mathbb{E}_{\boldsymbol{x}}\quantity[\left(y^{\star}(\boldsymbol{x})-\frac{\boldsymbol{x}\cdot\boldsymbol{w}}{\sqrt{d}}\right)^{2}]\,,$
(4)

where the average is over the sample distribution $P_{\rm
input}(\boldsymbol{x})$, and
$y^{\star}=\boldsymbol{x}\cdot\boldsymbol{w}^{\star}/\sqrt{d}$ is the
noiseless teacher-generated label. Thus, in this setting, cross-validation on
clean datasets is equivalent to cross-validation on the mean-square error with
the ground truth. We will be interested in characterising the behaviour of the
generalisation error and the estimation error when both the dimension $d$ and
the number of training samples $n$ are large, $d,n\to\infty$ with finite
sample complexity $\alpha=n/d$.

##  IV Characterisation of the performance of BO and ERM estimators

##### Bayes Optimal estimator.

The BO estimator is defined as the estimator achieving, on average over the
data distribution, the lowest value of a given performance metric. Define the
posterior distribution as

$P(\boldsymbol{w}\mid\mathcal{D})=\frac{1}{\mathcal{Z}_{\mathcal{D}}}P(\boldsymbol{w})P(\mathcal{D}\mid\boldsymbol{w})\,,$
(5)

where $\mathcal{Z}_{\mathcal{D}}$ is a normalisation factor, $\mathcal{D}$ is
the training dataset, and in our case
$P(\boldsymbol{w})=\mathcal{N}(0,\mathbbm{1}_{d})$. It is well known [27] that
the optimal estimator error, as defined in eq. (3), is achieved from the
posterior average

$\hat{\boldsymbol{w}}^{\rm BO}_{\rm
estim}(\mathcal{D})=\mathbb{E}_{\boldsymbol{w}|\mathcal{D}}\quantity[\boldsymbol{w}]$
(6)

Instead, the optimal estimator for the generalisation error, as defined in
eq.(2) is given by (see Appendix A.1)

$\hat{\boldsymbol{w}}^{\rm BO}_{\rm
gen}(\mathcal{D})=\mathbb{E}_{\boldsymbol{w}|\mathcal{D}}\quantity[\mathbb{E}_{\bar{\boldsymbol{x}},\bar{y}|\boldsymbol{w}}\quantity[\bar{y}\,\bar{\boldsymbol{x}}]]\,,$
(7)

where $\mathbb{E}_{\boldsymbol{w}|\mathcal{D}}\quantity[\cdot]$ is the
posterior average, and
$\mathbb{E}_{\bar{\boldsymbol{x}},\bar{y}|\boldsymbol{w}}\quantity[\cdot]$ is
the average over a sample-label pair $(\bar{\boldsymbol{x}},\bar{y})$
generated with the data model eq. (1) using $\boldsymbol{w}$ as a teacher
vector.

In the high-dimensional limit, one can characterise analytically the behaviour
of both the BO generalisation error and the BO estimation error. This was
first rigorously proven in [26], whose results we adapt to our situation to
characterise the information theoretic Bayesian bound. In [26] the authors
showed that in generalized linear models, the value of the estimation error
concentrates on

$E_{\mathrm{estim}}^{\mathrm{BO}}=1-q_{\mathrm{b}}\,,$ (8)

while the value of the generalisation error concentrates on (here we specify
their equations to the particular model of data used in the present paper)

$E_{\mathrm{gen}}^{\mathrm{BO}}=1+\epsilon(\beta^{2}-1)-q_{\mathrm{b}}(1+\epsilon(\beta-1))^{2}+\Delta_{\text{eff}}\,,$
(9)

where
$\Delta_{\text{eff}}=(1-\epsilon)\Delta_{\text{\tiny{IN}}}+\epsilon\Delta_{\text{\tiny{OUT}}}$
and the order parameter $q_{\mathrm{b}}$ — which can be interpreted as the
normalised overlap between the posterior average and the teacher vector — can
be found by solving the system of two equations

$\displaystyle q_{\mathrm{b}}$
$\displaystyle=\frac{\hat{q}_{\mathrm{b}}}{1+\hat{q}_{\mathrm{b}}}\,,\quad\hat{q}_{\mathrm{b}}=\alpha\int\mathbb{E}_{\xi}\left[\mathcal{Z}_{\mathrm{out}}\left(y,q_{\mathrm{b}}^{1/2}\xi,\rho_{\mathrm{w}^{\star}}-q_{\mathrm{b}}\right)f_{\mathrm{out}}\left(y,q_{\mathrm{b}}^{1/2}\xi,\rho_{\mathrm{w}^{\star}}-q_{\mathrm{b}}\right)^{2}\right]\,dy,$
(10)

with $f_{\mathrm{out}}$ and $\mathcal{Z}_{\mathrm{out}}$ defined as

$\begin{split}\mathcal{Z}_{\mathrm{out}}(y,\omega,V)&=\mathbb{E}_{z\sim\mathcal{N}(\omega,V)}\left[P_{\mathrm{out}}\quantity(y|z)\right]\,,\quad
f_{\mathrm{out}}(y,\omega,V)=\partial_{\omega}\log(\mathcal{Z}_{\mathrm{out}}(y,\omega,V))\,.\end{split}$
(11)

As $\alpha\to\infty$, $q_{b}\to 1$, showing that for large sample complexity
the BO estimation goes to zero, while the BO generalisation error does not, as
a non-zero excess error remains. For this reason, in the following we will
often work with the excess generalisation error, defined as

$E_{\mathrm{gen}}^{\mathrm{excess}}(\boldsymbol{w})=E_{\mathrm{gen}}(\boldsymbol{w})-E_{\mathrm{gen}}^{\mathrm{BO}}(\alpha=+\infty)=E_{\mathrm{gen}}(\boldsymbol{w})-\epsilon(1-\epsilon)(1-\beta)^{2}-\Delta_{\text{eff}}\,.$
(12)

Notice that, in practical scenarios, computing the posterior average, or more
generally providing an algorithm that efficiently reaches Bayes optimal
performance — for example, in polynomial time in $n,d$ — is challenging. In
this specific case, an efficient algorithm to compute the BO estimator
accurately in the considered limit on large $n,d$ and fixed ratio $n/d$
exists. It is the generalised approximate message passing (GAMP) algorithm
[28, 29, 26], which we will use to plot numerical simulations for BO errors in
the following (see Appendix B.3 for details).

We will use the BO estimator as a baseline to compare the performance of ERM
with the selected losses to the optimal one in the presence of outliers.
Specifically, eq. (8) and eq. (9) provide the information-theoretical lower
bounds for the estimation error and generalisation error, respectively.

##### Empirical risk minimisation.

Define the risk function

$\mathcal{R}(\boldsymbol{w})=\sum_{i=1}^{n}\mathcal{L}(y_{i},\hat{y}_{i}(\boldsymbol{w}))+\lambda\>r(\boldsymbol{w})\,,$
(13)

where
$\hat{y}_{i}(\boldsymbol{w})=\boldsymbol{x}_{i}\cdot\boldsymbol{w}/\sqrt{d}$
are the labels predicted by the model on each training sample, $y_{i}$ are the
noisy labels of the training dataset, $\mathcal{L}$ is a convex function of
$\hat{y}$, $r$ a convex function of $\boldsymbol{w}$ and $\lambda\geq 0$.
Minimising the risk function $\mathcal{R}(\boldsymbol{w})$ provides an
estimate $\hat{\boldsymbol{w}}$ of the teacher vector
$\boldsymbol{w}^{\star}$, and this estimate is uniquely defined thanks to
convexity of the risk.

The risk function is a sum of a data-dependent term $\mathcal{L}(y,\hat{y})$
and of a regularisation function $r(\boldsymbol{w})$. For the data-dependent
loss function $\mathcal{L}$, we will consider three convex losses: $\ell_{2}$,
$\ell_{1}$ and Huber loss. The $\ell_{2}$ loss
$\ell_{2}(y,\hat{y})=\frac{1}{2}(y-\hat{y})^{2}$ is the most common choice in
regression and will provide a good baseline for more outlier-robust losses.
The $\ell_{1}$ loss $\ell_{1}(y,\hat{y})=\absolutevalue{y-\hat{y}}$ is the
simplest choice for outlier-robust regression, as it penalises less poor
interpolation of outliers w.r.t. the $\ell_{2}$ loss. The Huber loss

${\ell^{\rm
Huber}_{a}}(y,y^{\prime})=\begin{cases}\frac{1}{2}(y-y^{\prime})^{2}&\text{if
}\absolutevalue{y-y^{\prime}}<a\\\
a\absolutevalue{y-y^{\prime}}-\frac{a^{2}}{2}&\text{if
}\absolutevalue{y-y^{\prime}}\geq a\\\ \end{cases}$ (14)

has been proposed as a middle-ground between $\ell_{2}$ and $\ell_{1}$ [3]. It
is more resilient to outliers than the $\ell_{2}$ loss, while remaining smooth
and convex as the $\ell_{2}$ loss. Proper tuning of the scale parameter $a$
effectively produces an $\ell_{2}$ loss landscape for the non-outlying data
points and an $\ell_{1}$-like one for the outliers. For $a\to+\infty$ the
Huber loss reduces to the $\ell_{2}$ loss, while for $a\to 0$ there is a
straightforward mapping to a $\ell_{1}$ loss (one needs to rescale the total
loss by a factor $a$ to take this limit, obtaining an effective regularisation
$\lambda a$). While our analysis could be carried on for any choice of convex
regulariser $r(\boldsymbol{w})$, we will only study the $\ell_{2}$
regularisation $r(\boldsymbol{w})=\frac{1}{2}\|\boldsymbol{w}\|_{2}^{2}$.
Notice that due to our choices of scaling, for $\lambda$ to have an effect on
the norm of the learned weights in the limit $\alpha\to\infty$, one needs the
scaling $\lambda=\mathcal{O}(\alpha)$.

We will mostly consider optimally-tuned risk minimisation, meaning that we
will set $\lambda$, and $a$ in the case of the Huber loss, such as to minimise
either the generalisation or the estimation error. In practice, the former is
equivalent to cross-validation on a dataset statistically identical to the
training set, while the latter to cross-validation on a clean dataset, without
any noise or outliers.

We will perform numerical simulations for the errors of ERM estimators using,
for the $\ell_{2}$ loss, the explicit solution for its minimiser. For Huber
loss, we use the L-BFGS method [30] initialised on a random Gaussian weight.
For the $\ell_{1}$ loss we use its mapping onto Huber with small scale
parameter, setting $a=10^{-3}$.

Our first main result is the analytic characterisation of the generalisation
and estimation errors of the ERM estimator.

######  Theorem IV.1.

For the ERM estimator of the risk function eq. (13) with $\ell_{2}$
regularisation and $\lambda\geq 0$, under the data model defined in eq. (1)
with $\beta=0$ and in the high dimensional limit $n,d\to\infty$ with
$n/d=\alpha$ fixed, we have that the excess generalisation error and the
estimation error concentrate on

$\displaystyle E_{\mathrm{gen}}^{\mathrm{excess}}$
$\displaystyle=1-\epsilon+q-2m(1-\epsilon)\,,$ (15) $\displaystyle
E_{\mathrm{estim}}$ $\displaystyle=1+q-2m\,,$

where the values of $q$ and $m$ are the solutions of a system of six self-
consistent equations for the unknowns
$(m,q,\Sigma,\hat{m},\hat{q},\hat{\Sigma})$. Three equations are loss-
independent

$m=\frac{\hat{m}}{\lambda+\hat{\Sigma}}\,,\quad
q=\frac{\hat{m}^{2}+\hat{q}}{(\lambda+\hat{\Sigma})^{2}}\,,\quad\Sigma=\frac{1}{\lambda+\hat{\Sigma}}\,,$
(16)

while the other three are loss-dependent. For the $\ell_{2}$ loss we have that

$\displaystyle\hat{m}$
$\displaystyle=\frac{\alpha}{1+\Sigma}(1-\epsilon)\,,\quad\hat{\Sigma}=\frac{\alpha}{1+\Sigma}\,,$
(17) $\displaystyle\hat{q}$
$\displaystyle=\frac{\alpha}{(1+\Sigma)^{2}}\quantity[1+q+(1-\epsilon)\Delta_{\text{\tiny{IN}}}+\epsilon\Delta_{\text{\tiny{OUT}}}-\epsilon-2m\quantity(1-\epsilon)]\,.$

For the $\ell_{1}$ loss and for the Huber loss we have that:

$\displaystyle\hat{m}$
$\displaystyle=\frac{\alpha}{\nu}\quantity[(1-\epsilon)\erf\quantity(\chi_{\text{\tiny{IN}}})]\,,\quad\hat{\Sigma}=\frac{\alpha}{\nu}\quantity[(1-\epsilon)\erf\quantity(\chi_{\text{\tiny{IN}}})+\epsilon\erf\quantity(\chi_{\text{\tiny{OUT}}})]\,,$
(18) $\displaystyle\hat{q}$
$\displaystyle=\frac{\alpha}{\nu^{2}}\Bigg{[}(1-\epsilon)\left(\zeta_{\text{\tiny{IN}}}-\mu^{2}\right)\erf\quantity(\chi_{\text{\tiny{IN}}})+\epsilon\left(\zeta_{\text{\tiny{OUT}}}-\mu^{2}\right)\erf\quantity(\chi_{\text{\tiny{OUT}}})$
$\displaystyle\quad-\mu\sqrt{\frac{2}{\pi}}\left((1-\epsilon)\sqrt{\zeta_{\text{\tiny{IN}}}}\,e^{-\chi_{\text{\tiny{IN}}}^{2}}\right.+\left.\epsilon\sqrt{\zeta_{\text{\tiny{OUT}}}}\,e^{-\chi_{\text{\tiny{OUT}}}^{2}}\right)+\mu^{2}\Bigg{]}\,,$

where $\zeta_{\text{\tiny{IN}}}=\Delta_{\text{\tiny{IN}}}-2m+q+1$,
$\zeta_{\text{\tiny{OUT}}}=\Delta_{\text{\tiny{OUT}}}+q$, while the other
constants depend on the loss. For $\ell_{1}$ loss

$\mu^{\ell_{1}}=\nu^{\ell_{1}}=\Sigma\,,\quad\chi_{\text{\tiny{IN}}}^{\ell_{1}}=\Sigma/\sqrt{2\zeta_{\text{\tiny{IN}}}}\,,\quad\chi_{\text{\tiny{OUT}}}^{\ell_{1}}=\Sigma/\sqrt{2\zeta_{\text{\tiny{OUT}}}}\,,$
(19)

and for Huber loss

$\nu^{\text{Huber}}=\Sigma+1\,,\quad\mu^{\text{Huber}}=a\nu^{\text{Huber}}\,,\quad\chi_{\text{\tiny{IN}}}^{\text{Huber}}=\mu^{\text{Huber}}/\sqrt{2\zeta_{\text{\tiny{IN}}}}\,,\quad\chi_{\text{\tiny{OUT}}}^{\text{Huber}}=\mu^{\text{Huber}}/\sqrt{2\zeta_{\text{\tiny{OUT}}}}\,.$
(20)

The previous Theorem follows from Theorem 1 in [25], through a reduction to
our case and a computation of loss-specific quantities (namely their proximal
operators) that is provided in detail in Appendix A.2. We also provide a
similar result for $\beta>0$. The parameters $m$ and $q$ appearing in the
previous statement have a simple interpretation as the values around which the
teacher-student and student-student overlaps concentrate in high dimension

$m=\frac{1}{d}\boldsymbol{w}^{\star}\cdot\hat{\boldsymbol{w}}\,,\quad
q=\frac{1}{d}\norm{\hat{\boldsymbol{w}}}_{2}^{2}\,.$ (21)

In the case of Ridge regression we have that the equations can be solved
explicitly and get to an explicit form for the generalisation error. We show
this in Appendix A.2.1.

##  V Consistency (and lack thereof) of the generalisation & estimation error

We now study, for both generalisation and estimation error and at large sample
complexity $\alpha$, whether optimally-regularised ERM estimators achieve BO
performance or not, i.e. whether ERM is a consistent estimator.

We start by remarking that the large sample complexity asymptotics of the BO
generalisation error and of the BO estimation error is

$E_{\mathrm{estim}}^{\mathrm{BO}}=0+\mathcal{O}\left(\alpha^{-1}\right)\,,\quad
E_{\mathrm{gen}}^{\mathrm{BO}}=\epsilon(1-\epsilon)(1-\beta)^{2}+\Delta_{\rm
eff}+\mathcal{O}\left(\alpha^{-1}\right)$ (22)

which can be found by expanding around
$q_{b}=1+c/\alpha+\mathcal{O}\left(\alpha^{-2}\right)$ in eq. (10). We prove
the (fast) rate $1/\alpha$ for the errors in Appendix A.3. This rate is a
typical feature of BO estimators, see for example [24, 31].

Figure 1:  Generalisation and estimation errors as a function of the sample
complexity $\alpha$ for $\beta=0$, in two regimes where ERM do not always
achieve BO performance for large $\alpha$. (Left) Here we plot the excess
generalisation error for $\epsilon=0.6$, $\Delta_{\text{\tiny{IN}}}=1$ and
$\Delta_{\text{\tiny{OUT}}}=0.5$, and for all losses we use the optimal value
of $\lambda$. We plot two versions of the Huber loss, one with optimal scale
parameter $a$ and the other with fixed $a=1$. We see that only the error of
the $\ell_{2}$ and Huber loss with optimal scale $a$ (superimposed in the
plot) vanishes at large sample complexity, reaching BO performance. For fixed-
scale Huber loss, and $\ell_{1}$ loss, the error converges to a finite value.
The parameters used for the simulations (dots) of BO are $d=1000$ averaged
over $50$ samples, while for the ERM $d=200$ averaged over $1000$ samples.
(Right) Here we plot the estimation error for $\epsilon=0.3$,
$\Delta_{\text{\tiny{IN}}}=1.0$ and $\Delta_{\text{\tiny{OUT}}}=5.0$, and for
all losses we use optimal $(\lambda,a)$. We see that all losses converge to a
finite value of the estimation error, while the BO error goes to zero for
large sample complexity. The parameters used for BO are $d=4000$ averaged over
$10$ samples, while for the ERM $d=200$ aveaged over $1000$ samples.

We then consider ERM. We find that ERM always estimates correctly the
direction of the teacher weights. One could then expect that by adjusting the
$\ell_{2}$ regularisation parameter, ERM estimators could also achieve the BO
norm for the metric under consideration (generalisation/estimation error). We
show that this is not always possible if we restrict $\lambda\geq 0$ — a very
natural constraint given that for $\lambda<0$ the risk function for $\ell_{1}$
and Huber losses becomes non-convex and unbounded from below. This is due to
the fact that for some values of the parameters characterising the outliers’
distribution, the norm of ERM estimators is shorter than the BO one, requiring
a priori a negative regularisation to push it to larger norms.

######  Theorem V.1.

Under the same conditions of Theorem IV.1, as $\alpha\to\infty$ (taken after
the $n,d\to\infty$ limit), we have that

  * •

the normalised scalar product between ERM estimators (for the three losses
considered) with the teacher weights converges to one, regardless of the value
of the regulariser $\lambda$ and of the scale parameter $a$ of the Huber loss.

  * •

$\ell_{2}$ and optimally-scaled Huber ERM estimators achieve consistency with
respect to the generalisation error for all fixed values of $\lambda>0$ with
the BO rate $1/\alpha$ for all values of the parameters of the outliers’
distribution.

  * •

Optimally-regularised $\ell_{1}$ and suboptimally-scaled Huber ERM estimators
achieve consistency with respect to the generalisation error with the BO rate
$1/\alpha$ if and only if
$\Delta_{\text{\tiny{OUT}}}-\Delta_{\text{\tiny{IN}}}\geq(1-\beta)^{2}(2\epsilon-1)$
for $0<\beta<1$, and
$\Delta_{\text{\tiny{OUT}}}-\Delta_{\text{\tiny{IN}}}\leq(1-\beta)^{2}(2\epsilon-1)$
for $\beta>1$.

  * •

Optimally-regularised ERM estimators (for the three losses considered) achieve
consistency with respect to the estimation error with the BO rate $1/\alpha$
if and only if $\beta>1$.

We prove Theorem V.1 in Appendix A.4 by expanding at large $\alpha$ the self-
consistent equations presented in Theorem IV.1.

Thus, we proved that there exists some regions of parameters in which, for
both generalisation and estimation error, $\ell_{2}$ regularisation is not
sufficient to reach BO performance. In those regimes, consistency for the
generalisation error can always be recovered by cross-validating on a test set
to find the best norm, given that the angle between teacher and student
vanishes. Due to homogeneity properties of the three losses considered, the
rescaling can also be performed before training on the labels, as discussed in
Appendix A.4.6. The estimation error can be calibrated similarly, but we
stress that such a calibration requires either having access to an oracle
knowledge of the BO norm, or to a clean, unpolluted dataset to cross-validate
for the optimal norm. In the absence of such an ideal hold-out dataset, the
Huber and $\ell_{1}$ methods do not lead to consistent estimates of the
objective functions, and robust estimation is biased.

In Figure1 we show, for both the generalisation and estimation error, the
dependency of the error on the sample complexity $\alpha$ in two regimes were
some of the losses do not achieve consistency.

It is worth noticing that consistency could be achieved by carefully allowing
for negative $\ell_{2}$ regularisation. In the case of the $\ell_{2}$ loss, a
negative regularisation does not immediately break the convexity of the risk
function. Indeed, the risk function is convex up to
$\lambda=-(1-\sqrt{\alpha})^{2}$ (see Appendix A.5). Thus, for $\ell_{2}$ loss
consistency for the estimation error can be achieved without using an oracle
norm up to (in the large sample complexity limit) $\beta>1-\epsilon^{-1}$, by
allowing for negative lambda and retaining convexity of the risk function. On
the other hand, for $\ell_{1}$ and Huber loss, even an arbitrarily small
negative regularisation implies a non-convex and unbounded-from-below risk
function. It is reasonable to think though that an optimisation procedure
restricted to small enough norms would still result in a well-defined ERM
procedure, allowing for negative regularisation also in these case. We leave a
more thorough exploration of this effect for future work. We also mention that
negative regularisation has been considered in a series of recent works [32,
33, 19].

##  VI Dependency of the generalisation on the parameters

Figure 2:  Excess generalisation error (top panels) as a function of the
percentage of outliers’ $\epsilon$ (left) and of the outliers’ variance
$\Delta_{\text{\tiny{OUT}}}$ (right), along with the respective optimal hyper-
parameters $(\lambda,a)$ (bottom panels). In both plots $\alpha=10$,
$\beta=0$, $\Delta_{\text{\tiny{IN}}}=1$, and for all numerical simulations
(dots) have $d=200$, and are averaged over $1000$ samples. (Left) In this plot
$\Delta_{\text{\tiny{OUT}}}=5$. We observe that for small $\epsilon$ both the
$\ell_{2}$ and the Huber loss converge to BO performance, while $\ell_{1}$
does not. (Right) In this plot $\epsilon=0.3$. For
$\Delta_{\text{\tiny{OUT}}}<\Delta_{\text{\tiny{IN}}}$ we notice that the
Huber scale parameter diverges, so the Huber and $\ell_{2}$ loss become
identical. We observe that the BO estimator has a non-monotonic behaviour as a
function of $\Delta_{\text{\tiny{OUT}}}$, while ERM estimators are monotonic.
Moreover, we see that the optimal scale parameter for the Huber loss diverges
at a finite value of $\Delta_{\text{\tiny{OUT}}}$, signalling that the Huber
loss reduces to an $\ell_{2}$ loss below that value of noise.

We now investigate the dependence of on the parameters characterising the
amount of corruption caused by outliers. As in the previous section, we plot
the excess generalisation error, and we always consider optimally tuned
hyperparameters $\lambda$ and $a$. Again, we discuss the phenomenology arising
for a particular value for the parameters of the model, but we stress that the
phenomenology is qualitatively the same for other values of the parameters,
see Appendix C for additional experiments.

##### Dependency of the generalisation error on the amount of outliers.

We investigate the influence of the percentage of outliers $\epsilon$ on the
generalisation error in Figure 2 left panel. We identify two qualitatively
different regimes where the behaviour is different: the regime with small
$\epsilon$ where outliers are rare and vice-versa the regime with $\epsilon$
close to 1, where the samples that would be considered outliers in the
$\epsilon\to 0$ limit dominate the dataset.

We observe that in the small $\epsilon$ regime, for all losses, the ERM excess
generalisation error increases as the percentage of outliers increases. In the
limit of small $\epsilon$, the Huber loss performs similarly to $\ell_{2}$,
which is known to be Bayes optimal at $\epsilon=0$ if optimally regularised
[27]. In contrast, the $\ell_{1}$ loss converges to a higher value for the
generalisation error not reaching BO performances. For all finite values
$\epsilon$, ERMs fails to reach BO performance for all three choices of the
loss we studied.

By expanding the explicit asymptotic form for the generalisation error of ERM
with $\ell_{2}$ loss presented in Corollary A.1, we are able to provide an
analytical characterisation of the small $\epsilon$ behaviour of $\ell_{2}$
ERM.

######  Theorem VI.1.

As $\epsilon\to 0$ (taken after the $n,d\to\infty$ limit), we have that the
generalisation error of optimally-regularised $\ell_{2}$ loss and the optimal
regularisation satisfy

$E_{\mathrm{gen}}^{\ell_{2}}=E_{\mathrm{gen}}^{\mathrm{BO}}+\kappa\epsilon+\mathcal{O}(\epsilon^{2})\,,\quad\lambda=\Delta_{\text{\tiny{IN}}}+\lambda_{1}\epsilon+\mathcal{O}(\epsilon^{2})$
(23)

where $E_{\mathrm{gen}}^{\mathrm{BO}}$ is the Bayes optimal generalisation
error, and $\kappa$ and $\lambda_{1}$ are functions of the model parameters
that we provide in Appendix A.6.

In the $\epsilon\gtrsim 0.5$ regime, where outliers are dominant, we observe a
decrease in the excess generalisation error, which eventually vanishes as
$\epsilon$ approaches one. In this limit the majority of the dataset is
composed of data generated from a $\beta$-rescaled teacher, so that ERM will
be biased to smaller ($0<\beta<1$) or larger ($\beta>1$) norms than if trained
on a dataset with no outliers. As discussed previously, the optimal $\ell_{2}$
regularisation will counterbalance as much as possible this effect, aiming at
the BO norm. In the extreme case where $\beta=0$, outlier samples provide no
meaningful information on the ground truth and the BO estimator equals the
null vector $\boldsymbol{w}=0$. This is compatible with what we observe in the
figure, where $\lambda_{\mathrm{opt}}$ diverges as $\epsilon$ approaches one.

##### Dependency of the generalisation error on the variance of outliers.

In Figure 2, right panel, we consider the dependency of the generalisation
error as a function of the variance of the outliers
$\Delta_{\text{\tiny{OUT}}}$. We start by noticing that all ERM estimators
have a generalisation error that grows monotonically with
$\Delta_{\text{\tiny{OUT}}}$, while the BO generalisation error has a non-
monotonic behaviour, with a local maximum. This behaviour is also present for
$\beta>0$. Moreover, we notice that for large outliers’ variance the
$\ell_{2}$ loss performs worse than the more resilient $\ell_{1}$ and Huber
losses, as expected.

##### Performance difference between Huber and $\ell_{2}$ loss

We observe an interesting behaviour in comparing the Huber loss and the
$\ell_{2}$ loss. Observing Figure 2, right panel, we observe that as
$\Delta_{\text{\tiny{OUT}}}$ decreases the scale parameter value $a_{\rm opt}$
of the Huber loss diverges to infinity at a finite value of
$\Delta_{\text{\tiny{OUT}}}$. When $a\to+\infty$ the Huber loss converges to
the $\ell_{2}$ loss, as confirmed by the fact that their predicted
generalisation errors coincides (Figure 2 right). It is remarkable that
optimally-regularised Huber loss is not always strictly better than the
simpler $\ell_{2}$ loss in a whole region of model parameters.

We explore further this phase transition by plotting the difference between
the generalisation error of the optimal $\ell_{2}$ and Huber losses as a
function of both $\epsilon$ and $\Delta_{\text{\tiny{OUT}}}$ in Figure 3. We
remark that the two plots in Figure 2 correspond respectively to horizontal
and vertical slices in this plot. Also, we see that the difference in
performance between the two losses is relevant just in the region of large
$\Delta_{\text{\tiny{OUT}}}$ and intermediate $\epsilon$. Moreover, focusing
on the region of small $\epsilon$, large $\Delta_{\text{\tiny{OUT}}}$ and
large $\alpha$, i.e. for really rare and noisy outliers, Figure 3 suggests
that the more samples one has (larger $\alpha$), the less convenient it
becomes to use the Huber loss, at fixed value of all other parameters.

Figure 3:  Difference between the generalisation error of optimally-
regularised $\ell_{2}$ and Huber loss losses as a function of $\epsilon$ and
$\Delta_{\text{\tiny{OUT}}}$ for $\Delta_{\text{\tiny{IN}}}=1$, $\beta=0$.
Each shaded region denotes where the two losses have a non-zero difference,
i.e. where it is better to use the Huber loss, for different values of
$\alpha$ as indicated by the labels. Focusing on the region of small
$\epsilon$, large $\Delta_{\text{\tiny{OUT}}}$ and large $\alpha$, i.e. for
really rare and noisy outliers, the plot suggests that the more samples one
has (larger $\alpha$), the less convenient it becomes to use the Huber loss,
at fixed value of all other parameters.

##  VII Concluding remarks

We studied robust ERM in a model of data with outliers, focusing on three
commonly used losses ($\ell_{2}$-regularised $\ell_{2}$, $\ell_{1}$, and Huber
loss). The percentage of outliers, variance, and typical norm can be varied,
allowing for a thorough exploration of the effect of outliers on the
regression task. We derived analytically exact asymptotics for the
generalisation and estimation errors of ERM estimators trained on this data
model in the high-dimensional limit, where the number of samples $n$ and the
dimension $d$ are large with a fixed finite ratio, and we compared them with
BO estimators.

For the generalisation error, we found that optimally regularised ERM is
asymptotically consistent in the large sample complexity limit at the
condition of performing a norm calibration. For the estimation error, we
showed that due to a norm calibration mismatch, the consistency of the
estimator requires an oracle estimate of the optimal norm or the presence of a
cross-validation set not corrupted by the outliers. We also observed an
interesting phase transition as a function of the outlier parameters. In a
full region of the parameter space $(\epsilon,\Delta_{\text{\tiny{OUT}}})$,
the optimal scale parameter of the Huber loss $a_{\text{out}}$ diverges, and
the error of the ERM with Huber and $\ell_{2}$ losses coincides.

Our work opens the way to the study of more complicated settings, for instance
random features regression [34] or deep learning in the lazy/neural tangent
kernel regime [35], where one could study effects that are not directly
explorable in the vanilla generalised linear regression setting, such as over-
parametrization [36], benign overfitting [37], implicit regularisation [38],
and their interplay with outliers.

## Acknowledgements

We thank Lenka Zdeborová for fruitful discussions and insightful ideas, and
Yatin Dandi and Borja Rodriguez Mateos for useful discussions in the early
stages of this project. We acknowledge funding from the ERC under the European
Union’s Horizon 2020 Research and Innovation Program Grant Agreement
714608-SMiLe, as well as by the Swiss National Science Foundation grants SNSF
OperaGOST $200021\\_200390$ and SNSF TMPFP2_210012.

## References

  * [1] Peter J. Huber.  Robust Statistics, pages 1248–1251.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2011.
  * [2] G. E. P. Box and G. C. Tiao.  A bayesian approach to some outlier problems.  Biometrika, 55(1):119–129, 1968.
  * [3] Peter J. Huber.  Robust Estimation of a Location Parameter.  The Annals of Mathematical Statistics, 35(1):73 – 101, 1964.
  * [4] Peter J. Huber.  A Robust Version of the Probability Ratio Test.  The Annals of Mathematical Statistics, 36(6):1753 – 1758, 1965\.
  * [5] F Hampel.  Contributions to the theory of robust estimation [ph. d. thesis].  University of California, Berkeley, Calif, USA, 1968.
  * [6] Frank R Hampel.  A general qualitative definition of robustness.  The Annals of mathematical statistics the official journal of the Institute of Mathematical Statistics, 42(6), 12 1971.
  * [7] Kush Bhatia, Prateek Jain, Parameswaran Kamalaruban, and Purushottam Kar.  Consistent robust regression.  In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
  * [8] Arun Sai Suggala, Kush Bhatia, Pradeep Ravikumar, and Prateek Jain.  Adaptive hard thresholding for near-optimal consistent robust regression.  CoRR, abs/1903.08192, 2019.
  * [9] Prasad Raghavendra and Morris Yau.  List decodable learning via sum of squares.  In Proceedings of the 2020 ACM-SIAM Symposium on Discrete Algorithms (SODA), pages 161–180, 2020.
  * [10] James Berger and L. Mark Berliner.  Robust bayes and empirical bayes analysis with $\epsilon$-contaminated priors.  The Annals of Statistics, 14(2):461–486, 1986.
  * [11] Pierre C Bellec.  Out-of-sample error estimate for robust m-estimators with convex penalty, 2020.
  * [12] Derek Bean, Peter J. Bickel, Noureddine El Karoui, and Bin Yu.  Optimal m-estimation in high-dimensional regression.  Proceedings of the National Academy of Sciences of the United States of America, 110(36):14563–14568, 2013.
  * [13] Efthymios Tsakonas, Joakim Jaldén, Nicholas D. Sidiropoulos, and Björn Ottersten.  Convergence of the huber regression m-estimate in the presence of dense outliers.  IEEE Signal Processing Letters, 21(10):1211–1214, 2014.
  * [14] Ka-Veng Yuen and Gilberto A. Ortiz.  Outlier detection and robust regression for correlated data.  Computer Methods in Applied Mechanics and Engineering, 313:632–646, 2017.
  * [15] Arnak S. Dalalyan and Philip Thompson.  Outlier-robust estimation of a sparse linear model using l1-penalized huber’s m-estimator.  In Neural Information Processing Systems, 2019.
  * [16] Philippe Gagnon, Alain Desgagné, and Mylène Bédard.  A New Bayesian Approach to Robustness Against Outliers in Linear Regression.  Bayesian Analysis, 15(2):389 – 414, 2020.
  * [17] Noureddine El Karoui, Derek Bean, Peter J. Bickel, Chinghway Lim, and Bin Yu.  On robust regression with high-dimensional predictors.  Proceedings of the National Academy of Sciences, 110(36):14557–14562, 2013.
  * [18] David Donoho and Andrea Montanari.  High dimensional robust m-estimation: asymptotic variance via approximate message passing.  Probability Theory and Related Fields, 166(3):935–969, Dec 2016\.
  * [19] Trevor Hastie, Andrea Montanari, Saharon Rosset, and Ryan J. Tibshirani.  Surprises in high-dimensional ridgeless least squares interpolation.  The Annals of Statistics, 50(2):949 – 986, 2022.
  * [20] Christos Thrampoulidis, Ehsan Abbasi, and Babak Hassibi.  Precise error analysis of regularized $m$ -estimators in high dimensions.  IEEE Transactions on Information Theory, 64(8):5592–5628, 2018\.
  * [21] Liu Liu, Tianyang Li, and Constantine Caramanis.  High dimensional robust $m$-estimation: Arbitrary corruption and heavy tails, 2019.
  * [22] Hanwen Huang.  Asymptotic risk and phase transition of $l_{1}$-penalized robust estimator.  The Annals of Statistics, 48(5):3090 – 3111, 2020.
  * [23] Pierre C. Bellec, Yiwei Shen, and Cun-Hui Zhang.  Asymptotic normality of robust M-estimators with convex penalty.  Electronic Journal of Statistics, 16(2):5591 – 5622, 2022.
  * [24] Benjamin Aubin, Florent Krzakala, Yue Lu, and Lenka Zdeborová.  Generalization error in high-dimensional perceptrons: Approaching bayes error with convex optimization.  In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 12199–12210. Curran Associates, Inc., 2020.
  * [25] Bruno Loureiro, Cédric Gerbelot, Hugo Cui, Sebastian Goldt, Florent Krzakala, Marc Mézard, and Lenka Zdeborová.  Learning curves of generic features maps for realistic datasets with a teacher-student model*.  Journal of Statistical Mechanics: Theory and Experiment, 2022(11):114001, nov 2022.
  * [26] 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, mar 2019.
  * [27] Thomas M Cover and Joy A Thomas.  Information theory and statistics.  Elements of Information Theory, 1(1):279–335, 1991.
  * [28] Yoshiyuki Kabashima and Shinsuke Uda.  A bp-based algorithm for performing bayesian inference in large perceptron-type networks.  In Shoham Ben-David, John Case, and Akira Maruoka, editors, Algorithmic Learning Theory, pages 479–493, Berlin, Heidelberg, 2004. Springer Berlin Heidelberg.
  * [29] Sundeep Rangan.  Generalized approximate message passing for estimation with random linear mixing.  In 2011 IEEE International Symposium on Information Theory Proceedings, pages 2168–2172. IEEE, 2011.
  * [30] Dong C. Liu and Jorge Nocedal.  On the limited memory bfgs method for large scale optimization.  Mathematical Programming, 45(1):503–528, Aug 1989.
  * [31] David Haussler, Michael Kearns, Manfred Opper, and Robert Schapire.  Estimating average-case learning curves using bayesian, statistical physics and vc dimension methods.  Advances in Neural Information Processing Systems, 4, 1991.
  * [32] Denny Wu and Ji Xu.  On the optimal weighted $\ell_{2}$ regularization in overparameterized linear regression.  In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 10112–10123. Curran Associates, Inc., 2020.
  * [33] Dmitry Kobak, Jonathan Lomond, and Benoit Sanchez.  The optimal ridge penalty for real-world high-dimensional data can be zero or negative due to the implicit ridge regularization.  J. Mach. Learn. Res., 21(1), jun 2022.
  * [34] Ali Rahimi and Benjamin Recht.  Random features for large-scale kernel machines.  Advances in neural information processing systems, 20, 2007.
  * [35] Lenaic Chizat, Edouard Oyallon, and Francis Bach.  On lazy training in differentiable programming.  Advances in neural information processing systems, 32, 2019.
  * [36] 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\.
  * [37] 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.
  * [38] Suriya Gunasekar, Blake E Woodworth, Srinadh Bhojanapalli, Behnam Neyshabur, and Nati Srebro.  Implicit regularization in matrix factorization.  Advances in Neural Information Processing Systems, 30, 2017.
  * [39] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research, 12:2825–2830, 2011.

##  Appendix A Mathematical details

###  A.1 Claim: Bayes optimal estimator for the generalisation error

Define the generalisation error as in eq. (2)

$E_{\mathrm{gen}}(\boldsymbol{w})=\mathbb{E}_{(\bar{\boldsymbol{x}},\bar{y}|\boldsymbol{w}^{\star})}\quantity[\quantity(\bar{y}-\frac{\bar{\boldsymbol{x}}\cdot\boldsymbol{w}}{\sqrt{d}})^{2}]\,,$
(24)

where the average is over a sample-label pair $(\bar{\boldsymbol{x}},\bar{y})$
generated with the same distribution as the training set, namely with the same
teacher vector $\boldsymbol{w}^{\star}$. The posterior-averaged generalisation
error is

$\mathcal{R}(\boldsymbol{w})=\mathbb{E}_{\boldsymbol{w}^{\star}|\mathcal{D}}\quantity[\mathbb{E}_{(\bar{\boldsymbol{x}},\bar{y})|\boldsymbol{w}^{\star}}\quantity[\quantity(\bar{y}-\frac{\bar{\boldsymbol{x}}\cdot\boldsymbol{w}}{\sqrt{d}})^{2}]]$
(25)

where $\mathcal{D}$ is the observed training dataset,
$\mathbb{E}_{\boldsymbol{w}^{\star}|\mathcal{D}}\quantity[\cdot]$ is the
posterior average and $(\bar{\boldsymbol{x}},\bar{y})$ is a sample label pair.
Differentiating and setting the derivative to zero (in order to minimise the
posterior-averaged generalisation error) gives that the minimiser
$\hat{\boldsymbol{w}}$, i.e. the BO estimator, satisfies

$\hat{\boldsymbol{w}}^{\rm BO}_{\rm
gen}(\mathcal{D})=\mathbb{E}_{\boldsymbol{w}|\mathcal{D}}\quantity[\mathbb{E}_{\bar{\boldsymbol{x}},\bar{y}|\boldsymbol{w}}\quantity[\bar{y}\,\bar{\boldsymbol{x}}]]\,,$
(26)

which for our specific choice of data model reduces to

$\hat{\boldsymbol{w}}^{\rm BO}_{\rm
gen}(\mathcal{D})=\quantity(1-\epsilon+\beta\epsilon)\mathbb{E}_{\boldsymbol{w}|\mathcal{D}}\quantity[\boldsymbol{w}]\,.$
(27)

###  A.2 Theorem 4.1: self-consistent equations for the errors of ERM

The more general proof of the form of the self-consistent equations can be
found in [25]. Our case corresponds to taking
$\boldsymbol{\varphi}_{s}(\boldsymbol{x})=\boldsymbol{\varphi}_{t}(\boldsymbol{x})=\boldsymbol{x}$
which is equivalent to the setting presented in [24], of which we will follow
the notation.

In the following we will refer to the loss function as
$\mathcal{L}(\cdot,\cdot)$ and to the regularisation function as $r(\cdot)$.
We then define the general quantities that will appear in the fixed point
equations, they are:

$\displaystyle\mathcal{Z}_{\mathrm{out}^{\star}}(y,\omega,V)$
$\displaystyle=\mathbb{E}_{z\sim\mathcal{N}(\omega,V)}\left[P_{\mathrm{out}}\quantity(y|z)\right]\,,$
(28) $\displaystyle f_{\mathrm{out}^{\star}}(y,\omega,V)$
$\displaystyle=\partial_{\omega}\log(\mathcal{Z}_{\mathrm{out}^{*}}(y,\omega,V))\,,$
$\displaystyle f_{\mathrm{out}}(y,\omega,V)$
$\displaystyle=(\mathcal{P}_{V}[\mathcal{L}(y,\cdot)](\omega)-\omega)/V\,,$
$\displaystyle\mathcal{Z}_{\mathrm{w}^{*}}(\gamma,\Lambda)$
$\displaystyle=\mathbb{E}_{\omega\sim
P_{\mathbf{w}^{*}}}\left[e^{-\frac{1}{2}\Lambda\omega^{2}+\gamma\omega}\right]\,,$
$\displaystyle f_{\mathrm{w}^{*}}(\gamma,\Lambda)$
$\displaystyle=\partial_{\gamma}\log\mathcal{Z}_{\mathrm{w}^{*}}(\gamma,\Lambda)\,,$
$\displaystyle f_{\mathrm{w}}(\gamma,\Lambda)$
$\displaystyle=\mathcal{P}_{\Lambda^{-1}}[r(.)]\left(\Lambda^{-1}\gamma\right),$

where $\omega$ in the definition of
$\mathcal{Z}_{\mathrm{w}^{*}}(\gamma,\Lambda)$ is sampled with the teacher
weight distribution and $\mathcal{P}_{V}[f](x)$ is the proximal operator
defined as follows

$\mathcal{P}_{V}[f(y,\cdot)](x)=\arg\min_{z}\left[f(y,z)+\frac{1}{2V}(z-x)^{2}\right]\,.$
(29)

The self-consistent equations in general form then read

$\displaystyle m$
$\displaystyle=\mathbb{E}_{\xi}\left[\mathcal{Z}_{\mathrm{w}^{\star}}(\sqrt{\hat{\eta}}\xi,\hat{\eta})f_{\mathrm{w}^{\star}}(\sqrt{\hat{\eta}}\xi,\hat{\eta})f_{\mathrm{w}}\left(\hat{q}^{1/2}\xi,\hat{\Sigma}\right)\right],$
(30) $\displaystyle q$
$\displaystyle=\mathbb{E}_{\xi}\left[\mathcal{Z}_{\mathrm{w}^{\star}}(\sqrt{\hat{\eta}}\xi,\hat{\eta})f_{\mathrm{w}}\left(\hat{q}^{1/2}\xi,\hat{\Sigma}\right)^{2}\right],$
$\displaystyle\Sigma$
$\displaystyle=\mathbb{E}_{\xi}\left[\mathcal{Z}_{\mathrm{w}^{\star}}(\sqrt{\hat{\eta}}\xi,\hat{\eta})\partial_{\gamma}f_{\mathrm{w}}\left(\hat{q}^{1/2}\xi,\hat{\Sigma}\right)\right]\,,$
$\displaystyle\hat{m}$
$\displaystyle=\alpha\int\mathbb{E}_{\xi}\left[\mathcal{Z}_{\mathrm{out}^{\star}}(y,\sqrt{\eta}\xi,(1-\eta))\cdot
f_{\mathrm{out}^{\star}}\left(y,\sqrt{\eta}\xi,(1-\eta)\right)f_{\mathrm{out}}\left(y,q^{1/2}\xi,\Sigma\right)\right]\,dy,$
$\displaystyle\hat{q}$
$\displaystyle=\alpha\int\mathbb{E}_{\xi}\left[\mathcal{Z}_{\mathrm{out}^{\star}}\left(y,\sqrt{\eta}\xi,(1-\eta)\right)f_{\mathrm{out}}\left(y,q^{1/2}\xi,\Sigma\right)^{2}\right]dy\,,$
$\displaystyle\hat{\Sigma}$
$\displaystyle=-\alpha\int\mathbb{E}_{\xi}\left[\mathcal{Z}_{\mathrm{out}^{\star}}\left(y,\sqrt{\eta}\xi,(1-\eta)\right)\partial_{\omega}f_{\mathrm{out}}\left(y,q^{1/2}\xi,\Sigma\right)\right]dy\,.$

where we have denoted with $\mathbb{E}_{\xi}$ the expectation value over a
standard Gaussian variable $\xi\sim\mathcal{N}(0,1)$, and we defined
$\eta=m^{2}/q$ and $\hat{\eta}=\hat{m}^{2}/\hat{q}$.

In the main text we fixed the regularisation function to the $\ell_{2}$
regularisation, i.e. $r(\mathbf{w})=\frac{1}{2}\norm{\mathbf{w}}^{2}_{2}$ and
the teacher weights to be extracted from the distribution
$\boldsymbol{w}^{\star}\sim\mathcal{N}(0,\mathbbm{1}_{d})$. The quantities
$\mathcal{Z}_{\mathrm{w}^{\star}}$ and $f_{\mathrm{w}^{\star}}$ can be written
as

$\mathcal{Z}_{\mathrm{w}^{*}}(\gamma,\Lambda)=\frac{e^{\frac{\gamma^{2}}{2(\Lambda+1)}}}{\sqrt{\Lambda+1}}\,,\quad
f_{\mathrm{w}^{*}}(\gamma,\Lambda)=\frac{\gamma}{\Lambda+1}\,.$ (31)

Also we have that

$\mathcal{Z}_{\mathrm{w}}^{\ell_{2},\lambda}(\gamma,\Lambda)=\frac{e^{\frac{\gamma^{2}}{2(\lambda+\Lambda)}}}{\sqrt{\lambda+\Lambda}}\,,\quad
f_{\mathrm{w}}^{\ell_{2},\lambda}(\gamma,\Lambda)=\frac{\gamma}{\lambda+\Lambda}\,.$
(32)

This allows to write the first three self consistent equations as:

$m=\frac{\hat{m}}{\lambda+\hat{\Sigma}}\,,\quad
q=\frac{\hat{m}^{2}+\hat{q}}{(\lambda+\hat{\Sigma})^{2}}\,,\quad\Sigma=\frac{1}{\lambda+\hat{\Sigma}}\,.$
(33)

As explained in [24] these first three equations only depend on the
regularisation term and they will be thus fixed during our analysis.

Thus, to compute the performance of ERM with a given loss, we just need to
compute $\mathcal{Z}_{\mathrm{out}^{\star}},f_{\mathrm{out}^{*}}$ and
$f_{\mathrm{out}}$, and solve the corresponding system of self-consistent
equations. For our noise model $P_{\rm out}$ we have that

$\mathcal{Z}_{\mathrm{out}^{\star}}(y,\omega,V)=\mathcal{N}_{y}(\omega,V+\Delta_{\text{\tiny{IN}}})+\mathcal{N}_{y}(\beta\omega,V+\beta^{2}\Delta_{\text{\tiny{OUT}}})\,,$
(34)

and $f_{\mathrm{out}^{\star}}$ is given from the definition in eq. (28).

In the following subsections we specialise these equations to the $\ell_{2}$,
$\ell_{1}$ and Huber losses by computing $f_{\mathrm{out}}$, which is the only
loss-dependent part in the equations.

#####  $\ell_{2}$ loss case.

In this case, the $f_{\mathrm{out}}$ in eq. (28) reads

$f_{\mathrm{out}}^{\ell_{2}}(y,\omega,V)=\frac{y-\omega}{1+V}\,.$ (35)

Then, the remaining integrals can be solved explicitly as they are all
expectation of Gaussian variables, giving

$\displaystyle m$ $\displaystyle=\frac{\hat{m}}{\lambda+\hat{\Sigma}}\,,\quad
q=\frac{\hat{m}^{2}+\hat{q}}{(\lambda+\hat{\Sigma})^{2}}\,,\quad\Sigma=\frac{1}{\lambda+\hat{\Sigma}}\,,$
(36) $\displaystyle\hat{m}$
$\displaystyle=\frac{\alpha}{1+\Sigma}(1+\epsilon(\beta-1))\,,\quad\hat{\Sigma}=\frac{\alpha}{1+\Sigma}\,,$
$\displaystyle\hat{q}$
$\displaystyle=\frac{\alpha}{(1+\Sigma)^{2}}\quantity[1+q+\Delta_{\text{eff}}+\epsilon\quantity(\beta^{2}-1)-2m\quantity(1+\epsilon(\beta-1))]\,.$

For this specific case we have that one can find an explicit form for the
order parameters which is presented in Appendix A.2.1.

##### Useful Integrals.

To perform the double integration over $(y,\xi)$ for the last three equations
in eq. (30) we will use a linear change of variables. The relevant change of
variables is $(u,v)=(y-\beta\sqrt{\eta}\xi,y-\sqrt{q}\xi)$. The associated
Jacobian is constant and equal to $\absolutevalue{\beta\sqrt{\eta}-\sqrt{q}}$.
With this change of variable, we obtain explicit expressions for the following
integrals, that will be used below. The first integral integral is

$\displaystyle\mathcal{I}_{1}$
$\displaystyle(\Delta,\beta,\kappa,\Lambda,a,m,q)=\int_{\mathbb{R}^{2}}\differential{\xi}\differential{y}e^{-\frac{1}{2}\xi^{2}-\frac{1}{2}\frac{\quantity(\beta\sqrt{\eta}\xi-y)^{2}}{\Delta+\beta^{2}(1-\eta)}}\quantity(\beta\sqrt{\eta}\xi-y)\begin{cases}a&\text{if
}y-\sqrt{q}\xi>\kappa\\\ \frac{y-\sqrt{q}\xi}{\Lambda}&\text{if
}\absolutevalue{y-\sqrt{q}\xi}<\kappa\\\ -a&\text{if }y-\sqrt{q}\xi<-\kappa\\\
\end{cases}$ (37)
$\displaystyle=\frac{2\pi}{\Lambda}\quantity(\quantity(1-\eta)\beta^{2}+\Delta)^{\frac{3}{2}}\quantity[\frac{\sqrt{2}(\kappa-a\Lambda)e^{-\frac{\kappa^{2}}{2\quantity(\beta^{2}+\Delta-2\beta
m+q)}}}{\sqrt{\pi(q-2m\beta+\beta^{2}+\Delta)}}-\erf\quantity(\kappa/\sqrt{2(\beta^{2}+\Delta-2\beta
m+q)})]\,.$

The second integral is:

$\displaystyle\mathcal{I}_{2}$
$\displaystyle(\Delta,\beta,\kappa,\Lambda,a,m,q)=\int_{\mathbb{R}^{2}}\differential{\xi}\differential{y}e^{-\frac{1}{2}\xi^{2}-\frac{1}{2}\frac{\quantity(\beta\sqrt{\eta}\xi-y)^{2}}{\Delta+\beta^{2}(1-\eta)}}\begin{cases}\quantity(\frac{y-\sqrt{q}\xi}{\Lambda})^{2}&\text{if
}\absolutevalue{y-\sqrt{q}\xi}<\kappa\\\ a^{2}&\text{otherwise}\\\
\end{cases}$ (38)
$\displaystyle=\frac{2\pi}{\Lambda^{2}}\sqrt{\quantity(1-\eta)\beta^{2}+\Delta}\Bigg{[}a^{2}-\kappa\sqrt{2(q-2m\beta+\beta^{2}+\Delta)/\pi}\,e^{-\frac{\kappa^{2}}{2\quantity(\beta^{2}+\Delta-2\beta
m+q)}}$
$\displaystyle\qquad+\quantity(q-2m\beta+\beta^{2}+\Delta-a^{2})\erf\quantity(\kappa/\sqrt{2(\beta^{2}+\Delta-2\beta
m+q)})\Bigg{]}\,.$

The third integral is

$\displaystyle\mathcal{I}_{3}$
$\displaystyle(\Delta,\beta,\kappa,\Lambda,a,m,q)=\int_{\mathbb{R}^{2}}\differential{\xi}\differential{y}e^{-\frac{1}{2}\xi^{2}-\frac{1}{2}\frac{\quantity(y-\beta\sqrt{\eta}\xi)^{2}}{\Delta+\beta^{2}(1-\eta)}}\mathbb{I}_{\absolutevalue{y-\sqrt{q}\xi}<\kappa}$
(39)
$\displaystyle=2\pi\sqrt{\quantity(1-\eta)\beta^{2}+\Delta}\erf\quantity(\kappa/\sqrt{2(\beta^{2}+\Delta-2\beta
m+q)})/\Lambda\,,$

Using these integrals we can write explicitly the self-consistent equations,
leading to the equations presented in Theorem IV.1.

#####  $\ell_{1}$ loss case.

For the $\ell_{1}$ loss, the $f_{\mathrm{out}}$ in eq. (28) reads

$f_{\mathrm{out}}^{\ell_{1}}(y,\omega,V)=\begin{cases}\frac{1}{V}\quantity(y-\omega)&\text{if
}|\omega-y|\leq V\\\ \operatorname{sign}(\omega-y)&\text{otherwise
}\end{cases}\,.$ (40)

With the integration technique explained above we have that the explicit form
of the self-consistent equations is

$\displaystyle\hat{m}$
$\displaystyle=\frac{\alpha}{\Sigma}\quantity[(1-\epsilon)\erf\quantity(\chi_{\text{\tiny{IN}}})+\beta\epsilon\erf\quantity(\chi_{\text{\tiny{OUT}}})]\,,$
(41) $\displaystyle\hat{q}$
$\displaystyle=\frac{\alpha}{\Sigma^{2}}\Bigg{[}(1-\epsilon)\left(\zeta_{\text{\tiny{IN}}}-\mu^{2}\right)\erf\quantity(\chi_{\text{\tiny{IN}}})+\epsilon\left(\zeta_{\text{\tiny{OUT}}}-\mu^{2}\right)\erf\quantity(\chi_{\text{\tiny{OUT}}})$
$\displaystyle\quad-\mu\sqrt{\frac{2}{\pi}}\left((1-\epsilon)\sqrt{\zeta_{\text{\tiny{IN}}}}\,e^{-\chi_{\text{\tiny{IN}}}^{2}}\right.+\left.\epsilon\sqrt{\zeta_{\text{\tiny{OUT}}}}\,e^{-\chi_{\text{\tiny{OUT}}}^{2}}\right)+\mu^{2}\Bigg{]}\,,$
$\displaystyle\hat{\Sigma}$
$\displaystyle=\frac{\alpha}{\Sigma}\quantity[(1-\epsilon)\erf\quantity(\chi_{\text{\tiny{IN}}})+\epsilon\erf\quantity(\chi_{\text{\tiny{OUT}}})]\,,$

where, as in the main text, we defined
$\zeta_{\text{\tiny{IN}}}=\Delta_{\text{\tiny{IN}}}-2m+q+1$,
$\zeta_{\text{\tiny{OUT}}}=\Delta_{\text{\tiny{OUT}}}+\beta^{2}+q-2\beta m$
and

$\mu^{\ell_{1}}=\nu^{\ell_{1}}=\Sigma\,,\quad\chi_{\text{\tiny{IN}}}^{\ell_{1}}=\Sigma/\sqrt{2\zeta_{\text{\tiny{IN}}}}\,,\quad\chi_{\text{\tiny{OUT}}}^{\ell_{1}}=\Sigma/\sqrt{2\zeta_{\text{\tiny{OUT}}}}\,,$
(42)

##### Huber loss case.

For the Huber loss, the $f_{\mathrm{out}}$ in eq. (28) reads

$f_{\mathrm{out}}^{\mathrm{Huber}}(y,\omega,V)=\begin{cases}-a&\text{if
}y-\omega>a(V+1)\\\ \frac{y-\omega}{V+1}&\text{if }|y-\omega|<a(V+1)\\\
a&\text{otherwise}\end{cases}\,.$ (43)

we have that the integrals in eq. (30) become

$\displaystyle\hat{m}$
$\displaystyle=\frac{\alpha}{\Sigma+1}\quantity[(1-\epsilon)\erf\quantity(\chi_{\text{\tiny{IN}}})+\beta\epsilon\erf\quantity(\chi_{\text{\tiny{OUT}}})]\,,$
(44) $\displaystyle\hat{q}$
$\displaystyle=\frac{\alpha}{\quantity(\Sigma+1)^{2}}\Bigg{[}(1-\epsilon)\left(\zeta_{\text{\tiny{IN}}}-\mu^{2}\right)\erf\quantity(\chi_{\text{\tiny{IN}}})+\epsilon\left(\zeta_{\text{\tiny{OUT}}}-\mu^{2}\right)\erf\quantity(\chi_{\text{\tiny{OUT}}})$
$\displaystyle\quad-\mu\sqrt{\frac{2}{\pi}}\left((1-\epsilon)\sqrt{\zeta_{\text{\tiny{IN}}}}\,e^{-\chi_{\text{\tiny{IN}}}^{2}}\right.+\left.\epsilon\sqrt{\zeta_{\text{\tiny{OUT}}}}\,e^{-\chi_{\text{\tiny{OUT}}}^{2}}\right)+\mu^{2}\Bigg{]}\,,$
$\displaystyle\hat{\Sigma}$
$\displaystyle=\frac{\alpha}{\Sigma+1}\quantity[(1-\epsilon)\erf\quantity(\chi_{\text{\tiny{IN}}})+\epsilon\erf\quantity(\chi_{\text{\tiny{OUT}}})]\,$

where $\zeta_{\text{\tiny{IN}}}=\Delta_{\text{\tiny{IN}}}-2m+q+1$,
$\zeta_{\text{\tiny{OUT}}}=\Delta_{\text{\tiny{OUT}}}+\beta^{2}+q-2\beta m$
and

$\nu^{\text{Huber}}=\Sigma+1\,,\quad\mu^{\text{Huber}}=a\nu^{\text{Huber}}\,,\quad\chi_{\text{\tiny{IN}}}^{\text{Huber}}=\mu^{\text{Huber}}/\sqrt{2\zeta_{\text{\tiny{IN}}}}\,,\quad\chi_{\text{\tiny{OUT}}}^{\text{Huber}}=\mu^{\text{Huber}}/\sqrt{2\zeta_{\text{\tiny{OUT}}}}\,.$
(45)

####  A.2.1 Corollary: explicit form of the error for $\ell_{2}$

######  Corollary A.1.

In the Ridge regression case, for a generic value of $\beta$, with fixed
regularisation $\lambda\geq 0$ the generalisation error has the explicit form

$\displaystyle E_{\mathrm{gen}}=$ $\displaystyle\
1+(\beta^{2}-1)\epsilon+q-2m(1+(\beta-1)\epsilon)\,,$ (46) $\displaystyle m=$
$\displaystyle\ \frac{2\alpha\Gamma}{p+t+1}\,,$ $\displaystyle q=$
$\displaystyle\
\frac{4\alpha\quantity[\alpha\Gamma^{2}(p+t-3)+\Lambda(p+t+1)]}{(p+t+1)\quantity[p^{2}-2c+t^{2}+2t(p+1)+1]}\,,$

where
$\Delta_{\text{eff}}=(1-\epsilon)\Delta_{\text{\tiny{IN}}}+\epsilon\Delta_{\text{\tiny{OUT}}}$,
$p=\alpha+\lambda$, $c=\alpha-\lambda$, $t=\sqrt{(p-1)^{2}+4\lambda}$,
$\Gamma=1+\epsilon(\beta-1)$ and
$\Lambda=1+\Delta_{\text{eff}}+\epsilon(\beta^{2}-1)$.

###### Proof.

The form for the excess generalisation error is the same as in Theorem IV.1,
so we look for a solution for the self-consistent equations, eqs. (36). The
equations can be simplified if one defines $\Gamma=1+\epsilon(\beta-1)$ and
$\Lambda=1+\Delta_{\text{eff}}+\epsilon\quantity(\beta^{2}-1)$ with
$\Delta_{\text{eff}}=(1-\epsilon)\Delta_{\text{\tiny{IN}}}+\epsilon\Delta_{\text{\tiny{OUT}}}$.
To solve them we first notice that the equations for $\Sigma$ and
$\hat{\Sigma}$ only depend on themselves and thus can be solved, getting:

$\Sigma=\frac{1-\alpha-\lambda+t}{2\lambda}\,,\quad\hat{\Sigma}=\frac{-1+\alpha-\lambda+t}{2}\,,$
(47)

where here and in the following we have that
$t=\sqrt{4\lambda+(-1+\alpha+\lambda)^{2}}$. Plugging these values in the
equations for $m$ and $\hat{m}$ gives

$m=\frac{2\alpha\Gamma}{\alpha+\lambda+t+1}\,,\quad\hat{m}=-\frac{2\alpha\Gamma\lambda}{-\alpha+\lambda+t+1}$
(48)

Finally, plugging everything in the equations for $q$ and $\hat{q}$, we obtain

$\displaystyle q=$
$\displaystyle\frac{4\alpha\left(\alpha\Gamma^{2}(p+t-3)+\Lambda(p+t+1)\right)}{(p+t+1)\left(-2c+(p+t)^{2}+2t+1\right)}$
(49) $\displaystyle\hat{q}=$
$\displaystyle\frac{4\alpha\lambda^{2}}{(2\lambda-p+t+1)^{2}}\left(\frac{4\alpha\left(\alpha\Gamma^{2}(p+t-3)+\Lambda(p+t+1)\right)}{(p+t+1)\left(-2c+(p+t)^{2}+2t+1\right)}+\Lambda-\frac{4\alpha\Gamma^{2}}{p+t+1}\right)$

where we defined $p=\alpha+\lambda$ and $c=\alpha-\lambda$. Upon rearranging
of terms we have the same form for $q$ as in the Corollary. ∎

###  A.3 The rate of BO is $1/\alpha$

The equations for the BO are presented in the main text eq. (10). If we
suppose that $q_{\mathrm{b}}=1-c/\alpha+\mathcal{O}(\alpha^{-2})$ and
$\hat{q}_{\mathrm{b}}=\hat{c}\alpha+\mathcal{O}(1)$ we have that

$q_{\mathrm{b}}=\frac{\hat{q}_{\mathrm{b}}}{1+\hat{q}_{\mathrm{b}}}=1-\frac{1}{\hat{c}\alpha}$
(50)

thus the relation between the expansion coefficients which is $\hat{c}=1/c$.
The coefficient $\hat{c}$ is found from the first order expansion of the self-
consistent equations and it is equal to

$\hat{c}=\int\differential{\xi}\differential{y}e^{-\frac{1}{2}\xi^{2}}\frac{\left(\frac{\epsilon(\beta\xi-y)e^{-\frac{(y-\beta\xi)^{2}}{2\Delta_{\text{\tiny{OUT}}}}}}{\sqrt{2\pi}\Delta_{\text{\tiny{OUT}}}}+\frac{(\epsilon-1)(y-\xi)e^{-\frac{(y-\xi)^{2}}{2\Delta_{\text{\tiny{IN}}}}}}{\sqrt{2\pi}\Delta_{\text{\tiny{IN}}}}\right)^{2}}{\frac{\epsilon
e^{-\frac{(y-\beta\xi)^{2}}{2\Delta_{\text{\tiny{OUT}}}}}}{\sqrt{2\pi}\sqrt{\Delta_{\text{\tiny{OUT}}}}}-\frac{(\epsilon-1)e^{-\frac{(y-\xi)^{2}}{2\Delta_{\text{\tiny{IN}}}}}}{\sqrt{2\pi}\sqrt{\Delta_{\text{\tiny{IN}}}}}}$
(51)

Note that even if we supposed a faster decay for the second term of
$q_{\mathrm{b}}$ we would find anyway that
$\hat{q}_{\mathrm{b}}\propto\alpha$.

###  A.4 Theorem 5.1: large sample size limit for the self-consistent
equations for ERM

In this Appendix we derive the asymptotic value of the generalisation and
estimation errors of ERM for $\alpha\to\infty$. We will first consider the
case for the $\ell_{2}$ loss where the explicit solution of the order
parameters gives a more direct approach to study the above mentioned limit.
After we will consider the $\ell_{1}$ and Huber loss with a framework
together.

#####  $\ell_{2}$ loss.

In this case we can leverage Corollary A.1 to get the large $\alpha$ limit of
the order parameters and thus of the generalisation error. We notice that the
solution presented in Corollary A.1 are formal solution that do not suppose
anything on the scalings in $\alpha$ of the other parameters, i.e. $\lambda$.
We will consider the case where $\lambda$ is constant in $\alpha$ and see that
this suffices to have BO rates for the generalisation error. From eq. (46) we
have by expanding

$q=\Gamma^{2}+\frac{\Lambda-\Gamma^{2}(2\lambda+1)}{\alpha}+\frac{\Gamma^{2}\left(3\lambda^{2}-1\right)-2\lambda\Lambda+\Lambda}{\alpha^{2}}+\mathcal{O}\quantity(\frac{1}{\alpha^{3}})\,,$
(52)

and

$m=\Gamma-\frac{\Gamma\lambda}{\alpha}+\frac{\Gamma(\lambda-1)\lambda}{\alpha^{2}}+\mathcal{O}\quantity(\frac{1}{\alpha^{3}})\,,$
(53)

with the same definition of $\Gamma=1+\epsilon(\beta-1)$ and
$\Lambda=1+\Delta_{\text{eff}}+\epsilon(\beta^{2}-1)$. By the knowledge of the
order parameters and the formula for the generalisation and estimation error
one has that

$\displaystyle
E_{\mathrm{gen}}=\frac{\Delta_{\text{eff}}-(\beta-1)^{2}(\epsilon-1)\epsilon}{\alpha}$
(54)
$\displaystyle+\frac{2(\beta-1)^{2}\lambda(\epsilon-1)\epsilon+((\beta-1)\lambda\epsilon+\lambda)^{2}-\left((\beta-1)^{2}(\epsilon-1)\epsilon\right)-2\Delta_{\text{eff}}\lambda+\Delta_{\text{eff}}}{\alpha^{2}}+\mathcal{O}(\alpha^{-3})\,,$
$\displaystyle
E_{\mathrm{estim}}=(\beta-1)^{2}\epsilon^{2}+\frac{(\beta-1)\epsilon((2\lambda+1)(\epsilon-1)-\beta(2\lambda\epsilon+\epsilon-1))+\Delta_{\text{eff}}}{\alpha}$
$\displaystyle+\frac{\lambda((\beta-1)\epsilon(\lambda(3(\beta-1)\epsilon+4)-2\beta)+\lambda)-\left((\beta-1)^{2}(\epsilon-1)\epsilon\right)-2\Delta_{\text{eff}}\lambda+\Delta_{\text{eff}}}{\alpha^{2}}+\mathcal{O}(\alpha^{-3})\,.$

Notice that in the estimation error the value of the plateu depends only on
the noise parameters $\epsilon$ and $\beta$ but not on the noise variances.

#####  $\ell_{1}$ and Huber loss

We make the following ansatz for the leading order form of the following
parameters in $\alpha$

$m=m_{0}\,,\qquad
q=q_{0}\,,\qquad\Sigma=\frac{\Sigma_{0}}{\alpha}\,,\qquad\hat{m}=\hat{m}_{0}\alpha\,,\qquad\hat{q}=\hat{q}_{0}\alpha\,,\qquad\hat{\Sigma}=\hat{\Sigma}_{0}\alpha\,,$
(55)

where $m_{0},q_{0},\Sigma_{0},\hat{m}_{0},\hat{q}_{0}$ and $\hat{\Sigma}_{0}$
are constants not depending on $\alpha$ in the limit. We use the following
ansatz for $\lambda$ as a function of $\alpha$

$\lambda=\lambda_{0}+\lambda_{1}\alpha\,,$ (56)

again with $\lambda_{1}$ and $\lambda_{2}$ finite in the large $\alpha$ limit.
We now show that this ansatz leads to a system of self-consistent equations
involving only $m_{0},q_{0},\Sigma_{0},\hat{m}_{0},\hat{q}_{0}$ and
$\lambda_{1},\lambda_{0}$ and additionaly find the conditions for which the
first order term in $\lambda$ vanishes. This system characterises at leading
order the behaviour of the errors of ERM.

From the loss-independent self-consistent equations presented in Theorem IV.1

$m=\frac{\hat{m}}{\lambda+\hat{\Sigma}}\,,\quad
q=\frac{\hat{m}^{2}+\hat{q}}{(\lambda+\hat{\Sigma})^{2}}\,,\quad\Sigma=\frac{1}{\lambda+\hat{\Sigma}}\,$
(57)

we obtain at large $\alpha$ the equations

$m_{0}=\frac{\hat{m}_{0}}{\lambda_{1}+\hat{\Sigma}_{0}}\,,\quad
q_{0}=\frac{\hat{m}_{0}^{2}}{(\lambda_{1}+\hat{\Sigma}_{0})^{2}}=m_{0}^{2}\,,\quad\Sigma_{0}=\frac{1}{\lambda_{1}+\hat{\Sigma}_{0}}\,$
(58)

Notice that the expressions are independent of $\lambda_{0}$, which means that
for $\lambda$ to have an effect on the generalisation error it needs to have a
non-vanishing coefficient for the term linear in $\alpha$, i.e.
$\lambda_{1}\neq 0$. At leading order, the loss-independent constants defined
in Theorem IV.1 read

$\zeta_{\text{\tiny{IN}},0}=\Delta_{\text{\tiny{IN}}}+(m_{0}-1)^{2}\,,\qquad\zeta_{\text{\tiny{OUT}},0}=\Delta_{\text{\tiny{OUT}}}+(m_{0}-\beta)^{2}\,.$
(59)

#####  $\ell_{1}$ loss.

For the $\ell_{1}$ loss we have at leading order

$\chi_{\text{\tiny{IN}}}=\chi_{\text{\tiny{OUT}}}=\mu=\nu=\mathcal{O}(1/\alpha)\,,$
(60)

giving for the three loss-dependent self-consistent equations

$\displaystyle\hat{m}_{0}$
$\displaystyle=\sqrt{\frac{2}{\pi}}\left(\frac{1-\epsilon}{\sqrt{\zeta_{\text{\tiny{IN}},0}}}+\frac{\beta\epsilon}{\sqrt{\zeta_{\text{\tiny{OUT}},0}}}\right)\,,$
(61) $\displaystyle\hat{\Sigma}_{0}$
$\displaystyle=\sqrt{\frac{2}{\pi}}\left(\frac{1-\epsilon}{\sqrt{\zeta_{\text{\tiny{IN}},0}}}+\frac{\epsilon}{\sqrt{\zeta_{\text{\tiny{OUT}},0}}}\right)\,,$
$\displaystyle\hat{q}_{0}$ $\displaystyle=1\,.$

##### Huber loss.

For the Huber loss we will consider the scale parameter $a$ to be $\alpha$
independent. Thus we have at leading order

$\chi_{\text{\tiny{IN}}}=a/\sqrt{2\zeta_{\text{\tiny{IN}},0}}\,,\quad\chi_{\text{\tiny{OUT}}}=a/\sqrt{2\zeta_{\text{\tiny{OUT}},0}}\,,\quad\mu=a\,,\quad\nu=1\,,$
(62)

giving for the three loss-dependent self-consistent equations

$\displaystyle\hat{m}_{0}$
$\displaystyle=(1-\epsilon)\erf\quantity(\frac{a}{\sqrt{2\zeta_{\text{\tiny{IN}},0}}})+\beta\epsilon\erf\quantity(\frac{a}{\sqrt{2\zeta_{\text{\tiny{OUT}},0}}})\,,$
(63) $\displaystyle\hat{\Sigma}_{0}$
$\displaystyle=(1-\epsilon)\erf\quantity(\frac{a}{\sqrt{2\zeta_{\text{\tiny{IN}},0}}})+\epsilon\erf\quantity(\frac{a}{\sqrt{2\zeta_{\text{\tiny{OUT}},0}}})\,,$
$\displaystyle\hat{q}_{0}$ $\displaystyle=\
(1-\epsilon)\left(\zeta_{\text{\tiny{IN}},0}-a^{2}\right)\erf\left(\frac{a}{\sqrt{2\zeta_{\text{\tiny{IN}},0}}}\right)+\epsilon\left(\zeta_{\text{\tiny{OUT}},0}-a^{2}\right)\erf\left(\frac{a}{\sqrt{2\zeta_{\text{\tiny{OUT}},0}}}\right)\,.$

At leading order, the generalisation and estimation error read

$\displaystyle E_{\mathrm{gen}}$
$\displaystyle=1+\epsilon(\beta^{2}-1)+m_{0}^{2}-2m_{0}(1+\epsilon(\beta-1))+\Delta_{\text{eff}}=$
(64)
$\displaystyle=(m_{0}-1)^{2}-2m_{0}\epsilon(\beta-1)+\epsilon(\beta^{2}-1)+\Delta_{\text{eff}}\,,$
$\displaystyle E_{\mathrm{estim}}$
$\displaystyle=1+q_{0}-2m_{0}=(m_{0}-1)^{2}\,.$

####  A.4.1 Rates of ERM

Before deriving any conclusion from the set of self-consistent equations that
we just derived in the large $\alpha$ limit, it is worth considering what is
the scaling of the first subleading term for $m$ and $q$, as it determines the
scaling of the first subleading term for the generalisation and estimation
errors.

It is easy to see that a general expansion in integer powers of alpha for the
order parameters (so more general of the leading order ansatz of eq. (55))
still satisfies the self-consistent equations, as the order parameters enter
the equations only through $\mathcal{C}^{\infty}$ functions. We just need to
check that our functions are indeed smooth. This is contingent on the fact
that $\Sigma_{0}\neq 0$ and $\hat{\Sigma}_{0}+\lambda_{1}\neq 0$. The first
condition can just assumed to be true because of eq. (58). Indeed, for
$\Sigma_{0}$ to be zero we need $\hat{\Sigma}_{0}$ to diverge, which can’t be
since $\zeta_{\text{\tiny{IN}},0}>0$, $\zeta_{\text{\tiny{OUT}},0}>0$ as long
as $\Delta_{\text{\tiny{IN}}}>0$, $\Delta_{\text{\tiny{OUT}}}>0$ (equation
(59)). For the second one we can simply notice that as long as $\epsilon\leq
1$, $\hat{\Sigma}_{0}>0$.

Thus, the subleading order of $m$ and $q$, and thus of the errors, is at least
of order $\mathcal{O}(\alpha^{-1})$.

On the other hand, at the values of the parameters for which we have that the
error of ERM converges to the BO performance, we have that the subleading
order of the errors, is at most of order $\mathcal{O}(\alpha^{-1})$, which is
the BO scaling.

Thus, whenever ERM is consistent, it also achieves BO rates. A more detailed
computation would be needed to asses whether ERM is identical to BO at order
$\mathcal{O}(\alpha^{-1})$, but we believe that this will not be the case in
general.

####  A.4.2 Angle between ERM minimiser and teacher goes to zero

The angle between the student and the teacher weights is given in concentrates
in high-dimension onto [24]

$\theta_{\boldsymbol{w}^{\star},\hat{\boldsymbol{w}}}=\frac{1}{\pi}\arccos\left(\frac{m}{\sqrt{q}}\right)$
(65)

but from the loss-independent equations eq. (58) we have $q\to m^{2}$, so that
$\theta\to 0$ for large $\alpha$.

####  A.4.3 Consistency of generalisation error

In order for the generalisation error to be consistent, we need to impose that
the excess generalisation error vanishes, i.e. we require that

$E_{\mathrm{gen}}=(m_{0}-1)^{2}-2m_{0}\epsilon(\beta-1)+\epsilon(\beta^{2}-1)+\Delta_{\text{eff}}\overset{!}{=}\epsilon(1-\epsilon)(1-\beta)^{2}+\Delta_{\text{eff}}$
(66)

Which results in $m_{0}=\Gamma=1+\epsilon(\beta-1)$. This imposes a condition
on $\lambda_{1}$ because of eq. (58) which reads

$\lambda_{1}=\frac{\hat{m}_{0}}{\Gamma}-\hat{\Sigma}_{0}$ (67)

This needs to be compatible with $\lambda\geq 0$, which at large $\alpha$
implies $\lambda_{1}\geq 0$. Firstly we will consider the $\ell_{2}$ loss. In
this case we have $\hat{m}_{0}=\Gamma$ and $\hat{\Sigma}_{0}=1$, so
$\lambda_{1}=0$ meaning that the value of $\lambda$ doesn’t need to diverge
with $\alpha$ to have consistency of $E_{\mathrm{gen}}$. This result is
compatible to the one obtained by looking at the explicit solution of the
order parameters as in Corollary A.1, where we have found that even with a
finite value of $\lambda$ in the limit $\alpha\to\infty$ the generalisation
error has optimal rate.

For the other two losses the inequality in eq. (67) is equivalent to

$(1-\beta)(1-\epsilon)\epsilon\left[\erf\quantity(\frac{a}{\sqrt{2\zeta_{\text{\tiny{IN}},0}}})-\erf\quantity(\frac{a}{\sqrt{2\zeta_{\text{\tiny{OUT}},0}}})\right]\geq
0$ (68)

For $\beta<1$, considering that $\epsilon\in(0,1)$, remembering the
definitions in eq. (59), this condition is equivalent to

$\Delta_{\text{\tiny{OUT}}}-\Delta_{\text{\tiny{IN}}}\geq(1-\beta)^{2}(2\epsilon-1)\,.$
(69)

If $\beta>1$ the condition in eq. (69) flips the order of the condition.

We remark that in the limit $a\to\infty$ the equations is always satisfied.
This tells us that whenever L2 is better then Huber, the optimisation over $a$
of Huber will lead to $a\to\infty$, i.e. the Huber loss becomes the $\ell_{2}$
loss.

####  A.4.4 Consistency of estimation error

In order for the estimation error to be consistent, we need to impose that the
estimation error vanishes for $\alpha\to\infty$. This is the condition

$E_{\mathrm{estim}}=(m_{0}-1)^{2}\overset{!}{=}0\,,$ (70)

which is equivalent to $m_{0}=1$. As in Appendix A.4.3 we use eq. (58) and
obtain the value of $\lambda_{1}$ as

$\lambda_{1}=(\beta-1)\epsilon\erf\quantity(\frac{a}{\sqrt{2\zeta_{\text{\tiny{OUT}},0}}})\,,$
(71)

for the Huber loss. The other two losses can be obtained by setting
$a\to\infty$ for the $\ell_{2}$, and $a\to 0$ with the appropriate rescaling
of $\lambda_{1}$ for $\ell_{1}$. Both limits do not alter the overall sign of
$\lambda_{1}$, which just depends on $\beta-1$. This implies that the
estimator is can only be consistent in estimation for $\beta>1$ since we have
the constraint $\lambda\geq 0$.

####  A.4.5 Consistency of estimation error for $\ell_{2}$ with negative
regularisation

In Appendix A.5 we show that the Ridge regression risk can endure negative
regularisation without losing strict convexity for all
$\lambda>-(1-\sqrt{\alpha})^{2}$.

Using eq. (71), in the limit $a\to\infty$ to recover the $\ell_{2}$ loss, we
see that the condition above for lambda implies, at leading order in $\alpha$
and for the $\ell_{2}$ loss that

$(\beta-1)\epsilon>-1\,.$ (72)

####  A.4.6 Oracle norm/label rescaling restores consistency

We notice that since the angle between the teacher and the student goes to
zero as shown in Appendix A.4.2, it suffices to rescale the student optimally
to achieve consistency. This can be achieved either by knowing the norm, or by
cross-validation on the appropriate validation set to find the best norm.

Notice that all losses satisfy some kind of homogeneity property, i.e.

${\ell^{\rm
Huber}_{a}}\quantity(c\,y,\frac{\boldsymbol{x}\cdot\boldsymbol{w}}{\sqrt{d}})=c^{2}\,{\ell^{\rm
Huber}_{c\,a}}\quantity(y,\frac{\boldsymbol{x}\cdot\boldsymbol{w}}{c\sqrt{d}})\,,$
(73)
$\ell_{1}\quantity(c\,y,\frac{\boldsymbol{x}\cdot\boldsymbol{w}}{\sqrt{d}})=c\,\ell_{1}\quantity(y,\frac{\boldsymbol{x}\cdot\boldsymbol{w}}{c\sqrt{d}})\,,$
(74)
$\ell_{2}\quantity(c\,y,\frac{\boldsymbol{x}\cdot\boldsymbol{w}}{\sqrt{d}})=c^{2}\,\ell_{2}\quantity(y,\frac{\boldsymbol{x}\cdot\boldsymbol{w}}{c\sqrt{d}})\,.$
(75)

For this reason, the rescaling of the ERM estimator can be produced
equivalently by rescaling the labels of the training set, and the
hyperparamters, appropriately.

###  A.5 Claim: Ridge regression risk is still convex for not-too-large
negative regularisation

Consider the Ridge regression risk

$\mathcal{R}_{\mathrm{ridge}}(\boldsymbol{w})=\frac{1}{2}\sum_{i=1}^{n}\norm{y_{i}-\frac{\boldsymbol{x}_{i}\cdot\boldsymbol{w}}{\sqrt{d}}}^{2}_{2}+\frac{\lambda}{2}\norm{\boldsymbol{w}}^{2}_{2}\,.$
(76)

This is a quadratic form in the $d$-dimensional vector $\boldsymbol{w}$, with
quadratic part

$\mathcal{R}_{\mathrm{ridge}}^{\mathrm{quad.}}(\boldsymbol{w})=\frac{1}{2}\sum_{a,b=1}^{d}{w}^{a}{w}^{b}\left(\frac{\alpha}{n}\sum_{i=1}^{n}{x}_{i}^{a}{x}_{i}^{b}+\lambda\delta_{ab}\right)\overset{\mathrm{def}}{=}\frac{1}{2}\sum_{a,b=1}^{d}{w}^{a}{w}^{b}H_{a,b}\,,$
(77)

where the Hessian of the quadratic form is the matrix
$H=X+\lambda\,\mathrm{Id}$. The matrix $X$ is a $d\times d$ Marchenko-Pastur
random matrix. Supposing $\alpha=n/d>1$, this matrix has minimal eigenvalue
concentrating on $(1-\sqrt{1/\alpha})^{2}$, giving minimal Hessian eigenvalue
equal to

$\text{minimal eigenvalue}=(1-\sqrt{\alpha})^{2}+\lambda\,.$ (78)

Thus, the Ridge regression risk is convex as long as the Hessian is positive
definite, meaning

$\forall\lambda\,:\,\lambda>-(1-\sqrt{\alpha})^{2}\,.$ (79)

###  A.6 Theorem 6.1: small outliers’ percentage $\epsilon$ limit of the error
of ERM with $\ell_{2}$ loss

We start from the knowledge of the explicit form of the order parameters in
the Ridge regression case, presented in Corollary A.1. If one expands $m$ and
$q$ up to second order in $\epsilon$ one has that:

$m=m_{0}+m_{1}\epsilon+m_{2}\epsilon^{2}+\mathcal{O}(\epsilon^{3})\,,\>q=q_{0}+q_{1}\epsilon+q_{2}\epsilon^{2}+\mathcal{O}(\epsilon^{3})\,,\>\lambda=\lambda_{0}+\lambda_{1}\epsilon+\lambda_{2}\epsilon^{2}+\mathcal{O}(\epsilon^{3})\,,$
(80)

with coefficients that are constants in $\epsilon$. The generalisation error
up to second order in $\epsilon$ is

$\displaystyle E_{\mathrm{gen}}$
$\displaystyle={E_{\mathrm{gen}}}_{0}+{E_{\mathrm{gen}}}_{1}\epsilon+{E_{\mathrm{gen}}}_{2}\epsilon^{2}+\mathcal{O}(\epsilon^{3})$
(81)
$\displaystyle=1+q+\left(\beta^{2}-1\right)\epsilon+\Delta_{\text{eff}}-2m((\beta-1)\epsilon+1)$
$\displaystyle=(1-2m_{0}+q_{0}+\Delta_{\text{\tiny{IN}}})+((\beta-1)(\beta-2m_{0}+1)+\Delta_{\text{\tiny{OUT}}}-\Delta_{\text{\tiny{IN}}}-2m_{1}+q_{1})\epsilon$
$\displaystyle+(1-2(\beta-1)m_{1}-2m_{2})\epsilon^{2}+\mathcal{O}\left(\epsilon^{3}\right)$

This shows that the generalisation error has the correct expansion in powers
of $\epsilon$. We now need to find the optimal regularisation $\lambda$, and
plug it back in in the generalisation error to complete the proof.

To start, we notice that at $\epsilon=0$ the BO generalisation error can be
found explicitly, as the output channel become a single Gaussian, giving

$\displaystyle{E_{\mathrm{gen}}^{\mathrm{BO}}}_{0}$
$\displaystyle=\frac{1}{2}\left(\sqrt{(\alpha+\Delta_{\text{\tiny{IN}}}+1)^{2}-4\alpha}-\alpha+\Delta_{\text{\tiny{IN}}}+1\right)$
(82) $\displaystyle{E_{\mathrm{gen}}^{\ell_{2}}}_{0}$
$\displaystyle=\frac{1}{2}\left(\frac{\lambda_{0}(\alpha+\Delta_{\text{\tiny{IN}}})+\alpha(\alpha+\Delta_{\text{\tiny{IN}}}-2)+\Delta_{\text{\tiny{IN}}}+\lambda_{0}+1}{\sqrt{(\alpha+\lambda_{0}-1)^{2}+4\lambda_{0}}}-\alpha+\Delta_{\text{\tiny{IN}}}+1\right)$

where in the second line we compare with the zero-th order term of $\ell_{2}$
ERM. The two generalisation errors match for
$\lambda_{0}=\Delta_{\text{\tiny{IN}}}$, giving us the zero-th order value of
the expansion of the optimal regularisation.

To find the next orders for the optimal regularisation $\lambda$ we compute
order-by-order the derivative of the generalisation error with respect to
$\lambda$. At first order we have

$\partialderivative{{E_{\mathrm{gen}}}_{1}}{\lambda_{1}}=\frac{2\alpha(\lambda_{0}-\Delta_{\text{\tiny{IN}}})}{\left(\alpha^{2}+2\alpha(\lambda_{0}-1)+(\lambda_{0}+1)^{2}\right)^{3/2}}=0$
(83)

which vanishes as soon as $\lambda_{0}=\Delta_{\text{\tiny{IN}}}$ as computed
previously. Looking at the derivative at second order in $\epsilon$ allows to
fix $\lambda_{1}$ as

$\lambda_{1}=\frac{\left(\alpha^{2}+2\alpha(\lambda_{0}-1)+(\lambda_{0}+1)^{2}\right)\left(\beta^{2}-2\beta(\lambda_{0}+1)-\Delta_{\text{\tiny{IN}}}+\Delta_{\text{\tiny{OUT}}}+2\lambda_{0}+1\right)}{\alpha^{2}+\alpha(3\Delta_{\text{\tiny{IN}}}-\lambda_{0}-2)+(\lambda_{0}+1)(3\Delta_{\text{\tiny{IN}}}-2\lambda_{0}+1)}\,.$
(84)

Thus the optimal order parameter can be written as
$\lambda_{\mathrm{opt}}=\Delta_{\text{\tiny{IN}}}+\lambda_{1}\epsilon+\mathcal{O}(\epsilon^{2})$
and the value of the error is

$\displaystyle E_{\mathrm{gen}}^{\ell_{2}}$
$\displaystyle={E_{\mathrm{gen}}^{\mathrm{BO}}}_{0}+\bigg{[}\frac{1}{2\kappa}\Big{(}2\alpha^{2}(\beta-1)+\alpha(\beta(\beta+2\Delta_{\text{\tiny{IN}}}-6)-3\Delta_{\text{\tiny{IN}}}+\Delta_{\text{\tiny{OUT}}}+5)$
(85)
$\displaystyle+(\Delta_{\text{\tiny{IN}}}+1)\left(\beta^{2}-\Delta_{\text{\tiny{IN}}}+\Delta_{\text{\tiny{OUT}}}-1\right)\Big{)}+\frac{1}{2}\left(-2\alpha(\beta-1)+\beta^{2}+\Delta_{\text{\tiny{IN}}}-\Delta_{\text{\tiny{OUT}}}-1\right)\bigg{]}\epsilon$
$\displaystyle+\mathcal{O}(\epsilon^{2})\,,$

with the definitions of eq. (82) and

$\displaystyle\kappa$
$\displaystyle=\sqrt{\alpha^{2}+2\alpha(\Delta_{\text{\tiny{IN}}}-1)+(\Delta_{\text{\tiny{IN}}}+1)^{2}}\,.$
(86)

##  Appendix B Simulations details

###  B.1 Solving the self-consistent equations for ERM and BO

The self-consistent equations from Proposition IV.1 are written in a way
amenable to be solved via fixed-point iteration. From a random initialisation,
we iterate through both the hat and non-hat variable equations until the
maximum absolute difference between the order parameters in two successive
iterations falls below a tolerance of $10^{-9}$.

To speed-up convergence we use a damping scheme, updating each order parameter
at iteration $i$, designated as $x_{i}$, using
$x_{i}:=x_{i}\mu+x_{i-1}(1-\mu)$, with $\mu$ as the damping parameter
($\mu\in[0.7,0.9]$ was found most effective).

The Bayes-optimal (BO) fixed-point equations involve integrals which we could
not express explicitly, necessitating a numerical integrator. The integration
region must be subdivided into smaller regions to improve accuracy. In the
iteration of the self-consistent equations, the above-mentioned damping trick
is employed.

Once convergence is achieved for fixed $\lambda$ and possibly $a$, hyper-
parameters are optimised using the Nelder-Mead optimiser, a gradient-free
numerical minimisation procedure.

###  B.2 Numerical simulations for ERM

We obtain numerical simulations for the Empirical Risk Minimisation (ERM) by
sampling synthetic datasets using eq. (1). The ERM problem is solved as
follows:

  * •

for $\ell_{2}$ ERM is performed using the closed form solution. If we define
the data matrix $\Phi\in\mathbb{R}^{n,d}$ as $\Phi_{ij}=x_{i}^{j}$, then the
$\ell_{2}$ ERM problem is solved by:

$\hat{\boldsymbol{w}}=\quantity(\frac{\Phi^{\top}\Phi}{d}+\lambda\,\mathrm{Id})^{-1}\frac{\Phi^{\top}\boldsymbol{y}}{\sqrt{d}}$
(87)

  * •

For the Huber loss, ERM is performed using the routine of the Scikit-learn
library [39], which uses a L-BFGS optimiser [30].

  * •

For the $\ell_{1}$ loss we use its mapping onto Huber with small scale
parameter, setting $a=10^{-3}$. To obtain the same result as the Huber
regression with regularisation parameter $\lambda$ oen has to consider the
$\ell_{1}$ regression with regularisation parameter $a\lambda$.

For the comparison with synthetic data in Figure 1 we first find the optimal
parameters $(\lambda,a)$ for the generalisation error from the explicit
solution of the fixed point equations. These values are then plugged into the
ERM procedures described above.

To evaluate the estimation error it is pretty straightforward from the
definition in eq. (3). To evaluate the excess generalisation error from the
definition of eq. (2), we generated $n_{\mathrm{t}}$ sample dataset from the
same teacher and, given an estimate for the weights $\hat{\boldsymbol{w}}$ we
evaluated the following

$\frac{1}{n_{\mathrm{t}}}\sum_{i=1}^{n_{\mathrm{t}}}\quantity(y_{i}-\frac{\boldsymbol{x}\cdot\hat{\boldsymbol{w}}}{\sqrt{d}})^{2}-\frac{1}{n_{\mathrm{t}}}\sum_{i=1}^{n_{\mathrm{t}}}\quantity(y_{i}-\quantity(1-\epsilon+\beta\epsilon)\frac{\boldsymbol{x}\cdot\boldsymbol{w}^{\star}}{\sqrt{d}})^{2}\,,$
(88)

where the rescaling $(1-\epsilon+\epsilon\beta)$ is the one with minimal
generalisation error.

${{f_{\mathrm{out}}}_{\mu}^{t}}$${\gamma_{i}^{t}}$${\hat{\mathrm{w}}^{t+1}_{i}}$${\omega_{\mu}^{t+1}}$${{f_{\mathrm{out}}}_{\mu}^{t+1}}$${\gamma_{i}^{t+1}}$${\Lambda_{\mu}^{t}}$${c^{t+1}_{i}}$${V_{\mu}^{t+1}}$${\Lambda_{i}^{t+1}}$
Figure 4: Sketch of the dependence of the various quantities in GAMP
algorithm.

###  B.3 Numerical simulations for BO: the Generalised Approximate Message
Passing (GAMP) algorithm

We summarise here the Generalised Approximate Message Passing (GAMP)
algorithm, firstly introduced in [29]. We will follow the notation in [24].
Given loss and regularisation functions $f_{\mathrm{out}}$ and
$f_{\mathrm{w}}$, computed as in Appendix (A.2), we can estimate the student
weight $\hat{\boldsymbol{w}}$ of the associated ERM problem by iterating the
following fixed point system

$\begin{aligned}
{f_{\mathrm{out}}}_{\mu}^{t}&=f_{\mathrm{out}}(y_{\mu},\omega_{\mu}^{t},V_{\mu}^{t})\\\
\hat{\boldsymbol{w}}^{t+1}_{i}&=f_{\mathrm{w}}(\gamma^{t}_{i},\Lambda^{t}_{i})\\\
c^{t+1}_{i}&=\partial_{\gamma}f_{\mathrm{w}}(\gamma^{t}_{i},\Lambda^{t}_{i})\\\
\end{aligned}\qquad\qquad\qquad\begin{aligned}
\Lambda_{i}^{t}&=-\frac{1}{d}\sum_{\mu=1}^{n}X_{\mu
i}^{2}\partial_{\omega}{f_{\mathrm{out}}}_{\mu}^{t}\\\
\gamma_{i}^{t}&=\frac{1}{\sqrt{d}}\sum_{\mu=1}^{n}X_{\mu
i}{f_{\mathrm{out}}}_{\mu}^{t}+\Lambda_{i}^{t}\hat{\boldsymbol{w}}^{t}_{i}\\\
V_{\mu}^{t}&=\frac{1}{d}\sum_{i=1}^{d}X_{\mu i}^{2}\hat{c}^{t}_{i}\\\
\omega_{\mu}^{t}&=\frac{1}{\sqrt{d}}\sum_{\mu=1}^{n}X_{\mu
i}\hat{\boldsymbol{w}}^{t}_{i}-V^{t}_{\mu}{f_{\mathrm{out}}}_{\mu}^{t-1}\\\
\end{aligned}$ (89)

It has been shown [26] that if $f_{\mathrm{w}}=f_{\mathrm{w}^{\star}}$ and
$f_{\mathrm{out}}=f_{\mathrm{out}^{\star}}$ then GAMP provides a Bayes Optimal
estimation $\hat{\boldsymbol{w}}=\boldsymbol{w}_{\rm{BO}}$ as per Section IV.
We sketch in Figure 4 the dependence of the parameters of GAMP during one
iteration step going from iteration time $t$ to iteration time $t+1$.

##  Appendix C Exploration of parameters’ phenomenology

In this appendix we complement the figures in the main text by offering a more
comprehensive exploration of the dependence of parameters $\alpha$,
$\epsilon$, $\Delta_{\text{\tiny{IN}}}$, $\Delta_{\text{\tiny{OUT}}}$ and
$\beta$.

In most of the figures

we plot both the error (generalisation or estimation) in the upper panel and
the value of the optimal hyper parameters in the bottom panel. Hyperparameters
are plotted with colors matching the associated upper panle.

Gray lines correspond to the optimally tuned value of the scale parameter $a$
for Huber.

In Figure 5, 6, 7, 8 and 9 we explore the dependence on $\alpha$:

  * •

in Figures 5 and 6 we plot the optinally-tuned excess generalisation error as
a function of $\alpha$ for two choices of $\beta$. In figure 5 we displayed
the dependence on $\Delta_{\text{\tiny{OUT}}}$. We notice that the
phenomenology is qualitatively the same for the two different values of
$\beta$ ($\beta=0$ on the left column and $\beta=0.2$ on the right column). As
we expect from Theorem V.1. for large enough values of
$\Delta_{\text{\tiny{OUT}}}$ we achieve consistency with all losses, in
particular with a diverging $\lambda_{\mathrm{opt}}$. On the other hand if the
value of $\Delta_{\text{\tiny{OUT}}}$ is smaller than the threshold in V.1 no
choice of $\lambda\geq 0$ restores consistency in the large $\alpha$ limit. In
Figure 6 we show also that the same phenomenology holds when changing
$\epsilon$. As in Figure 5 we see that for certain choices of the parameters,
at large enough $\epsilon$ the estimator is not consistent.

  * •

In Figures 7 and 8, we consider the optimally-tuned estimation error. In the
figure 7 we focus on $\beta<1$ and see that estimation cannot be found for
positive values of $\lambda$. In figure 8 we have instead $\beta>1$. Here
consistent estimation can be achieved for an optimally tuned $\lambda$, but we
stress that this requires a clean dataset.

  * •

In Figure 9 we plotted the teacher-student angle as per eq. (65). We see that
the figures are consistent with the claims of Section A.4.2 since even for
non-optimised values of the parameters, the angle between the teacher and
student vector is going to zero with $\alpha$.

Figure 10 and Figure 11 explore the dependency of the generalisation error on
$\epsilon$ and $\Delta_{\text{\tiny{OUT}}}$. We have qualitatively the same
behaviour as in Figure 2 of the main text. Interestingly, the excess
generalisation error is not zero in the limit $\epsilon\to 1$ for $\beta\neq
0$ for the parameters chosen in figures. The outlier regression problem in
this limit is equivalent to regression without outlier, but with variance
$\Delta_{\text{\tiny{OUT}}}$.

Figure 12 shows quantitatively how the Huber loss performs better than the
$\ell_{2}$ loss in the regression task at hand when optimally tuned. This
complements Figure 3 of the main text, where we simply delimited the region
where the two losses are identical. Interestingly, the advantage of Huber over
$\ell_{2}$ grows the further away one is from the borders of this region.

Figure 5:  Excess generalisation error as a function of sample complexity for
different value of the noise models’ parameters. The hyper parameters,
$\lambda$ and optionally $a$, for the ERM procedures are optimised to obtain
the best generalisation error.

Figure 6:  Excess generalisation error as a function of sample complexity for
different value of the noise models’ parameters. In all the plots shown we
fixed $\Delta_{\text{\tiny{IN}}}=1$. The hyper parameters, $\lambda$ and
optionally $a$, for the ERM procedures are optimised to obtain the best
generalisation error.

Figure 7:  Estimation error as a function of sample complexity for different
value of the noise models’ parameters for $\beta<1$. In all the plots shown we
fixed $\Delta_{\text{\tiny{IN}}}=1$. The hyper parameters, $\lambda$ and
optionally $a$, for the ERM procedures are optimized to obtain the best
generalisation error.

Figure 8:  Estimation error as a function of sample complexity for different
value of the noise models’ parameters for $\beta>1$. In all the plots shown we
fixed $\Delta_{\text{\tiny{IN}}}=1$. The hyper parameters, $\lambda$ and
optionally $a$, for the ERM procedures are optimized to obtain the best
generalisation error.

Figure 9:  Angle between the teacher and the student — eq. (65) — as a
function of the sample complexity $\alpha$. The figures are generated for a
fixed value of regularisation parameter $\lambda=10^{-3}$ and in the case of
the Huber loss for $a=1$. For each plot we fixed $\Delta_{\text{\tiny{IN}}}=1$
and the other parameters chosen are indicated on top of the plot.

Figure 10:  Excess generalisation error as a function of $\epsilon$. In all
the plots shown we fixed $\alpha=10$. The hyper parameters, $\lambda$ and
optionally $a$, for the ERM procedures are optimized to obtain the best
generalisation error.

Figure 11:  Excess generalisation error as a function of
$\Delta_{\text{\tiny{OUT}}}$. In all the plots shown we fixed $\alpha=10$. The
hyper parameters, $\lambda$ and optionally $a$, for the ERM procedures are
optimized to obtain the best generalisation error.

Figure 12:  Difference between the Huber and $\ell_{2}$ generalisation error
as a function of $\Delta_{\text{\tiny{OUT}}}$ and $\epsilon$. Here we fix
$\Delta_{\text{\tiny{IN}}}=1$. Hyper-parameters for both losses are optimally
tuned. In the dark blue region the difference is identically zero.