# On Regression in Extreme Regions

Nathan Huet
LTCI, Telecom Paris
Institut Polytechnique de Paris
&Stephan Clémençon
LTCI, Telecom Paris
Institut Polytechnique de Paris
&Anne Sabourin
MAP5-UMR 8145
Université Paris Cité

###### Abstract

In the classic regression problem, the value of a real-valued random variable
$Y$ is to be predicted based on the observation of a random vector $X$, taking
its values in $\mathbb{R}^{d}$ with $d\geq 1$ say. The statistical learning
problem consists in building a predictive function
$\hat{f}:\mathbb{R}^{d}\to\mathbb{R}$ based on independent copies of the pair
$(X,Y)$ so that $Y$ is approximated by $\hat{f}(X)$ with minimum error in the
mean-squared sense. Motivated by various applications, ranging from
environmental sciences to finance or insurance, special attention is paid here
to the case of extreme (i.e. very large) observations $X$. Because of their
rarity, they contribute in a negligible manner to the (empirical) error and
the predictive performance of empirical quadratic risk minimizers can be
consequently very poor in extreme regions. In this paper, we develop a general
framework for regression in the extremes. It is assumed that $X$’s conditional
distribution given $Y$ belongs to a non parametric class of heavy-tailed
probability distributions. It is then shown that an asymptotic notion of risk
can be tailored to summarize appropriately predictive performance in extreme
regions of the input space. It is also proved that minimization of an
empirical and non asymptotic version of this ’extreme risk’, based on a
fraction of the largest observations solely, yields regression functions with
good generalization capacity. In addition, numerical results providing strong
empirical evidence of the relevance of the approach proposed are displayed.

_Keywords_ Statistical Learning Theory $\cdot$ Extreme Value Theory $\cdot$
Nonparametric Regression $\cdot$ Conditional Regular Variation

##  1 Introduction

Regression is a predictive problem of crucial importance in statistical
learning, covering a wide variety of applications. In the standard setup,
$(X,Y)$ is a pair of random variables defined on the same probability space
$(\Omega,\;\mathcal{A},\;\mathbb{P})$ with distribution $P$, where $Y$ is a
square integrable real-valued r.v. (the output) and $X$ is a random vector
with marginal distribution $\rho$ taking its values in some measurable space
$\mathcal{X}$ modelling some input information hopefully useful to predict
$Y$. The predictive learning problem consists in building, from a training
dataset $\mathcal{D}_{n}=\\{(X_{1},Y_{1}),\;\ldots,\;(X_{n},Y_{n})\\}$
composed of $n\geq 1$ independent copies of $(X,Y)$, a mapping
$f:\mathcal{X}\rightarrow\mathbb{R}$ in order to compute a ’good’ prediction
$f(X)$ for $Y$, with the quadratic risk

$R_{P}(f)=\mathbb{E}\left[\left(Y-f(X)\right)^{2}\right]$ (1)

as close as possible to that of $f^{*}(X)=\mathbb{E}[Y\mid X]$, which
obviously minimizes (1) over the space $L_{2}(\rho)$ of square integrable
functions of $X$:

$R_{P}^{*}:=R_{P}(f^{*})=\min_{f\in L_{2}(\rho)}R_{P}(f).$

A natural strategy consists in solving the Empirical Risk Minimization problem
(ERM in abbreviated form) $\min_{f\in\mathcal{F}}R_{\hat{P}_{n}}(f)$, where
$\mathcal{F}\subset L_{2}(\rho)$ is a closed and convex class of functions
sufficiently rich to include a reasonable approximant of $f^{*}$ and
$\hat{P}_{n}$ is a statistical version of the unknown distribution $P$,
typically the raw empirical distribution $(1/n)\sum_{i\leq
n}\delta_{(X_{i},Y_{i})}$, denoting by $\delta_{a}$ the Dirac mass at any
point $a$. The performance of predictive functions $\hat{f}_{n}$ obtained this
way, by least-square regression, has been extensively investigated in the
statistical learning literature Györfi et al. (2002); Massart (2007). Under
the assumption that the tails of the random pairs $(f(X),Y)$ are subgaussian
and appropriate complexity conditions are satisfied by the class
$\mathcal{F}$, confidence upper bounds for the excess of quadratic risk
$R_{P}(\hat{f}_{n})-R_{P}^{*}=\mathbb{E}[(Y-\hat{f}_{n}(X))^{2}\mid\mathcal{D}_{n}]-R_{P}^{*}$
have been established in Lecué and Mendelson (2016) by means of concentration
inequalities for empirical process Boucheron et al. (2013). We now place
ourselves in the situation where $\mathcal{X}=[0,\;+\infty)^{d}$. Observations
are considered as extreme when their norm exceeds some (asymptotically) large
threshold $t_{n}>0$. The threshold $t_{n}$ depends on the number $n\geq 1$ of
observations, since ’large’ should be naturally understood as large w.r.t the
vast majority of data observed. Hence, extreme observations are rare by nature
and severely underrepresented in the training dataset $\mathcal{D}_{n}$ with
overwhelming probability. Consequently, the impact of prediction errors in
extreme regions of the input space on the global regression error of
$\hat{f}_{n}$ is generally negligible. Indeed, the law of total probability
yields

$R_{P}(f)=\mathbb{P}\\{\|X\|>t_{n}\\}\mathbb{E}\left[\left(Y-f(X)\right)^{2}\mid\|X\|>t_{n}\right]\\\
+\mathbb{P}\\{\|X\|\leq
t_{n}\\}\mathbb{E}\left[\left(Y-f(X)\right)^{2}\mid\|X\|\leq t_{n}\right].$
(2)

Because the order of magnitude of $\mathbb{P}\\{\|X\|>t_{n}\\}$ (as well as
that of its empirical version with high probability) is extremely small, there
is no guarantee that the standard ERM strategy produces a predictive function
$\hat{f}$ that is nearly optimal in the extreme region
$\\{x:\;\|x\|>t_{n}\\}$, i.e. is such that the conditional error
$\mathbb{E}[(Y-\hat{f}(X))^{2}\mid\|X\|>t_{n}]$ is nearly minimum. However,
accurate prediction in extreme regions turns out to be crucial in certain
practical (safety) applications, in environmental sciences, dietary risk
analysis or finance/insurance for instance. It is thus the purpose of the
subsequent analysis to investigate the problem of building a prediction
function that asymptotically minimizes the first term of the decomposition
(2), and thus the quantity referred to as the quadratic conditional risk and
given by

$R_{t}(f):=\mathbb{E}\left[\left(Y-f(X)\right)^{2}\mid\|X\|>t\right]=R_{P_{t}}(f),$
(3)

denoting by $P_{t}$ the conditional distribution of $(X,Y)$ given $\|X\|>t$.
In order to develop a framework showing that empirical quadratic conditional
risk minimization leads to predictive rules with good generalization
capacities in extreme regions, (nonparametric) assumptions related to the tail
behavior of the conditional distribution of $X$ given $Y$ are required.
Multivariate regular variation hypotheses are very flexible in the sense that
they correspond to a large nonparametric class of (heavy-tailed)
distributions. They are frequently used in applications where the impact of
extreme observations should be enhanced, or not neglected at the minimum. This
hypothesis has been used to formulate unsupervised learning problems in the
extremes, anomaly detection in Thomas et al. (2017) and dimensionality
reduction in Goix et al. (2017). It is also used in Jalalzai et al. (2018) in
order to develop a framework for binary classification in extreme regions:
precisely, it is assumed therein that the two class distributions belong to
the maximal domain of attraction of multivariate extreme value distributions
with the same tail/shape index. Here we propose to work under a nonparametric
assumption related to the joint distribution of $(X,Y)$, ensuring in
particular that $X$’s conditional distribution given $Y$ belongs to the
maximal domain of attraction of a multivariate extreme value distribution with
probability one and guaranteeing that the conditional quadratic risk (3) may
converge to a risk functional referred to as asymptotic quadratic conditional
risk and describing asymptotic predictive performances in the extremes. Under
mild assumptions, we prove that a predictive rule using the angular
information only i.e. of the form $f(X)=f_{\Theta}(X/||X||)$, where
$f_{\Theta}$ is a real-valued function defined on the intersection of the
hypersphere and the positive orthant
$\mathbb{S}=\\{x\in[0,+\infty)^{d}:\;||x||=1\\}$, just like the minimizer of
the asymptotic conditional quadratic risk, learned by minimizing an empirical
version of (12) based on a fraction $k/n$ of the training dataset
$\mathcal{D}_{n}$, those corresponding to the largest $||X_{i}||$’s, is nearly
optimal w.r.t. the asymptotic conditional quadratic risk. Precisely,
nonasymptotic bounds for its excess of asymptotic conditional quadratic risk
are established. Beyond these theoretical guarantees, the performance of
empirical risk minimization in the extremes is supported by various numerical
experiments that have been carried out.

The paper is organized as follows. In Section 2, key notions pertaining to
multivariate extreme value theory (MEVT) are briefly recalled for clarity’s
sake and the probability framework we consider for regression in extreme
regions is described at length. The approach we propose for regression in the
extremes is detailed in Section 3, together with the dedicated validity
framework we develop. Illustrative experimental results are displayed in
Section 4, while some concluding remarks are collected in Section 5. Due to
space limitations, certain technical details are deferred to the Supplementary
Material.

##  2 Background and Preliminaries

We start with recalling some concepts in heavy-tail analysis, involved in the
formulation of the statistical problem given next. Here and throughout,
$(X,Y)$ is a pair of random variables defined on a probability space
$(\Omega,\;\mathcal{A},\;\mathbb{P})$ with distribution $P$, where $Y$ is
real-valued with marginal distribution $G$ and
$X=(X^{(1)},\;\ldots,\;X^{(d)})$ takes its values in $[0,\;+\infty)^{d}$,
$d\geq 1$, with marginal distribution $\rho$. We denote by $F_{j}$ the
(univariate) cumulative distribution function (cdf) of the component
$X^{(j)}$, $j=1,\;\ldots,\ d$, by $\mathds{1}\\{\mathcal{E}\\}$ is meant the
indicator function of any event $\mathcal{E}$, the integer part of any
$u\in\mathbb{R}$ is denoted by $\lfloor u\rfloor$ and $\|\cdot\|$ refers to a
norm on $\mathbb{R}^{d}$. Let
$E=[0,+\infty]^{d}\setminus\\{(0,\;\ldots,\;0)\\}$ be the punctured positive
orthant in $\mathbb{R}^{d}$ and denote by $\mathcal{B}(E)$ its Borel
$\sigma$-algebra. The boundary and the closure of any subset $B$ of the
topological space $\mathbb{R}^{d}$ are respectively denoted by $\partial B$
and $\bar{B}$ and we set $tB=\\{tx:\;x\in B\\}$ for all $t\in\mathbb{R}$. We
also consider the set $\mathbb{S}=\\{x\in[0,+\infty)^{d}:\;\|x\|=1\\}$ and
denote by $\mathcal{B}(\mathbb{S})$ its Borel $\sigma$-algebra.

###  2.1 Heavy Tails - Multivariate Regular Variation

The goal of heavy-tail analysis is to study phenomena that are not ruled by
averaging effects but determined by extreme values. To investigate the
behavior of a r.v. $X$ far from the center of its mass, an usual assumption is
that $X$’s distribution is multivariate regularly varying with tail index
$\alpha>0$, i.e. there exist a non zero Borel measure $\mu$ on $E$, finite on
all Borel measurable subsets of $E$ bounded away from zero and a regularly
varying111Recall that a positive Borel measurable function
$b:\mathbb{R}\to\mathbb{R}_{+}$ is said to be regularly varying with index
$\zeta\in\mathbb{R}$ iff $b(tu)/b(u)\to t^{\zeta}$ as $u\to+\infty$ for all
$t>0$. function $b(t)$ with index $\alpha$ such that

$b(t)\mathbb{P}\left\\{t^{-1}X\in B\right\\}\rightarrow\mu(B)\text{ as
}t\rightarrow+\infty,$ (4)

for any Borel measurable set $B\subset E$ bounded away from zero and s.t.
$\mu(\partial B)=0$. The limit measure $\mu$ is referred to as the exponent
measure of the r.v. $X$ and is provably homogeneous of degree $-\alpha$:
$\mu(tB)=t^{-\alpha}\mu(B)$ for all $t>0$ and Borel set $B\subset E$ bounded
away from the origin. One may refer to Resnick (2013) for alternative
formulations/characterizations of the regular variation property and its
application to MEVT. It follows from the homogeneity property that the
pushforward measure of $\mu$ by the polar coordinates transformation $x\in
E\mapsto(||x||,\Theta(x))$, where $\Theta(x):=x/||x||$ for all $x\in E$, is
the tensor product given by:

$\mu\left(\left\\{x\in E:\;||x||\geq r,\;\Theta(x)\in
A\right\\}\right)=r^{-\alpha}\Phi(A),$

for all $A\in\mathcal{B}(\mathbb{S})$ and $r\geq 1$, where $\Phi$ is a finite
positive measure on $\mathbb{S}$, referred to as the angular measure of the
heavy-tailed r.v. $X$. The regular variation assumption of the r.v. $X$
implies that the conditional distribution of $(\|X\|/t,\Theta(X))$ given
$\|X\|>t$ converges as $t\to+\infty$: for all $r\geq 1$ and
$A\in\mathcal{B}(\mathbb{S})$ with $\Phi(\partial A)=0$

$\mathbb{P}\Big{\\{}t^{-1}\|X\|\geq r,\Theta(X)\in
A\mid\|X\|>t\Big{\\}}\underset{t\rightarrow+\infty}{\longrightarrow}cr^{-\alpha}\Phi(A),$
(5)

where $c=\Phi(\mathbb{S})^{-1}=\mu(E\setminus\mathbb{B})^{-1}$ with
$\mathbb{B}:=\\{x\in E,\|x\|\leq 1\\}$. Hence, the radial and angular
components of the r.v. $X$ are asymptotically independent with standard Pareto
of parameter $\alpha$ and normalized angular measure $c\Phi$ as respective
asymptotic marginal distributions. The angular measure $\Phi$ describes
exhaustively the dependence structure of its components $X^{(j)}$ in the
extremes, i.e. the directions $\Theta(X)$ in which extremes occur with largest
probability.

Heavy-tailed models have been the subject of much attention in the machine-
learning literature. Among many other works, the regular variation assumption
is used in Ohannessian and Dahleh (2012) for rare event probability
estimation, in Bubeck et al. (2013), Carpentier and Valko (2014) or Achab et
al. (2017) in the context of stochastic bandit problems, in Goix et al. (2015)
for the statistical recovery of the dependence structure in the extremes, in
Goix et al. (2017) for dimensionality reduction in extreme regions and in
Brownlees et al. (2015) for predictive problems with heavy-tailed losses.

Here and throughout, for simplicity, we assume that $\alpha=1$ rather than
performing a standardization of $X$’s components to unit-Pareto margins, i.e.
replacing $X^{(j)}$ by $1/(1-F_{j}(X^{(j)}))$ for $j\in\\{1,\;\ldots,\ d\\}$,
see Remark 1 below.

###### Remark 1

(On marginal standardization) We point out that, by means of the marginal
standardization described above, the case of non-standard regular variation,
_e.g._ allowing for different tail indices $\alpha_{j}\neq 1$ for the
different components $X^{(j)}$, can also be reduced to the framework
considered here. However, the marginal transformations are unknown in
practice, just like the $F_{j}$’s and must be replaced with their empirical
counterparts, see Section 3. As will be discussed therein, this leads to
technical difficulties inherent to the dependence of the resulting pseudo-
standardized observations Clémençon et al. (2021).

In order to formulate the regression problem in the extremes, a specific,
though general/nonparametric, hypothesis related to the distribution of the
pair $(X,Y)$ is required, extending the usual multivariate regular variation
property. It is detailed and discussed in the following subsection.

###  2.2 Regression in the Extremes - The Framework

We now describe rigorously the framework we consider for regression in extreme
regions. For simplicity, we suppose that $Y$ is bounded through this paper.
This assumption can be naturally relaxed at the price of additional
technicalities.

###### Assumption 1

The r.v. $Y$ is bounded: there exists $M<+\infty$ s.t. $|Y|\leq M$ almost-
surely.

The hypothesis below, related to the asymptotic of $(X,Y)$’s conditional joint
distribution (suitably renormalized) given $\|X\|>t$ as $t\rightarrow+\infty$,
can be viewed as an extension of (4).

###### Assumption 2

There exist a transition kernel from $[-M,M]$ to $E$ given by
$(y,B)\in[-M,M]\times\mathcal{B}(E)\mapsto\mu_{y}(B)$, with $\mu_{y}$ a non-
null Radon measure on $\mathcal{B}(E)$ for every $y\in[-M,M]$, and a regularly
varying function $b$ with index $1$ such that

$\lim_{t\rightarrow+\infty}b(t)\mathbb{P}\left\\{t^{-1}X\in B\mid
Y=y\right\\}=\mu_{y}(B),$ (6)

for all $y\in[-M,M]$ and $B\in\mathcal{B}(E)$ bounded away from zero and s.t.
$\mu_{y}(\partial B)=0$.

In addition, $\sup_{t>1}b(t)\mathbb{P}\left\\{t^{-1}X\in B\mid Y\right\\}$ is
an integrable random variable.

For all $y\in[-M,M]$, the limit measure $\mu_{y}$ is referred to as the
conditional exponent measure of $X$ given $Y=y$. As can be shown by means of a
straightforward dominated convergence argument, under Assumption 2, for all
$B\in\mathcal{B}(E)$ bounded away from zero s.t. $\mu_{y}(\partial B)=0$, the
mapping $y\mapsto\mu_{y}(B)$ is integrable w.r.t. $G$ and $X$’s marginal
distribution is regularly varying with exponent measure
$\mu(A)=\lim_{t\rightarrow+\infty}b(t)\mathbb{P}\left\\{X\in
tA\right\\}=\int_{-M}^{M}\mu_{y}(A)dG(y)$. Alternative formulations of the
conditional regular variation hypothesis above are given in the Supplementary
Material. One may also consider the conditional angular measure of $X$ given
$Y=y$, denoted by $\Phi_{y}$, given by

$\frac{\mathbb{P}\\{\Theta(X)\in A,\|X\|\geq tr\mid
Y=y\\}}{\mathbb{P}\\{\|X\|\geq t\mid
Y=y\\}}\underset{t\rightarrow+\infty}{\longrightarrow}c_{y}r^{-1}\Phi_{y}(A),$
(7)

with $c_{y}=\Phi_{y}(\mathbb{S})^{-1}=\mu_{y}(E\setminus\mathbb{B})^{-1}$, for
all $A\in\mathcal{B}(\mathbb{S})$, $r\geq 1$ and $y\in[-M,M]$. In particular,
the two conditional limit measures are linked through the relation

$\mu_{y}(\\{x\in E,\|x\|\geq r,\Theta(x)\in A\\})=r^{-1}\Phi_{y}(A)$ (8)

for all $y\in[-M,M],r\geq 1$ and $A\in\mathcal{B}(\mathbb{S})$. Observe that

$\lim_{t\rightarrow\infty}\frac{\mathbb{P}\left\\{\Theta(X)\in A,\|X\|\geq
t\mid Y=y\right\\}}{\mathbb{P}\left\\{\|X\|>t\mid
Y=y\right\\}}=c_{y}\Phi_{y}(A),$

for all $y\in[-M,M]$, so that $c_{y}\Phi_{y}$ is the asymptotic conditional
probability distribution of $\Theta(X)$ given $Y=y$ and $\|X\|>t$ as
$t\to+\infty$. Notice also that $X$’s angular measure writes
$\Phi(A)=\int_{-M}^{M}\Phi_{y}(A)dG(y)$. Denote by
$G_{t}(y)=\mathbb{P}\\{Y\leq y\mid\|X\|>t\\}$ the conditional cdf of $Y$ given
$\|X\|>t$. As $t\to+\infty$, under Assumption 2, a dominated convergence
argument shows that it converges everywhere to the cdf

$G_{\infty}(y):=\Phi(\mathbb{S})^{-1}\int_{-M}^{y}\Phi_{u}(\mathbb{S})dG(u),\;\;y\in[-M,M].$
(9)

Let $P_{\infty}(dx\times dy)$ be the probability distribution on
$\mathcal{B}(E\setminus\mathbb{B})\times[-M,M]$ determined by its second
marginal distribution $G_{\infty}$ and $c_{y}\mu_{y}$ its conditional
distribution given $y$. Then, using (9), we have: $\forall
B\in\mathcal{B}(E\setminus\mathbb{B})$, $\forall y\in[-M,M]$,

$P_{\infty}\left\\{B\times[-M,y]\right\\}=\int_{-M}^{y}c_{u}\mu_{u}(B)dG_{\infty}(u)=\int_{-M}^{y}\frac{\mu_{u}(B)}{\Phi_{u}(\mathbb{S})}\frac{\Phi_{u}(\mathbb{S})}{\Phi(\mathbb{S})}dG(u)=\lim_{t\rightarrow+\infty}\mathbb{P}\Big{\\{}t^{-1}X\in
B,Y\leq y\mid\|X\|\geq t\Big{\\}}.$ (10)

The limit result above describes the asymptotic behavior of the conditional
distribution of $(X/t,Y)$ given $\|X\|>t$ as $t\rightarrow+\infty$. Here and
throughout, $(X_{\infty},Y_{\infty})$ shall denote a pair of r.v. with
distribution $P_{\infty}$. Observe that, under Assumptions 1 and 2, the r.v.
$Y_{\infty}$ is almost-surely bounded in amplitude by $M<+\infty$.

The example below describes a typical situation, where Assumption 1 and
Assumption 2 are satisfied. Refer also to the Supplementary Material for an
additional generic example (i.e. multiplicative noise model).

###### Example 1

(Additive noise model with heavy-tailed random design) Suppose that $X$ is a
heavy-tailed random vector, independent from a real-valued r.v.
$\varepsilon_{0}$, bounded and centered, modelling some noise, drawn from a
distribution with continuous density function $g_{0}$ w.r.t. the Lebesgue
measure in which the r.v.

$Y=f^{*}(X)+\varepsilon_{0}$ (11)

is observed, where $f^{*}:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is a Borel-
measurable and bounded mapping. Assume also that there exists a continuous
(necessarily bounded) function
$f^{*}_{\Theta}:\mathbb{S}\rightarrow\mathbb{R}$ s.t.

$\sup_{x\in E:\|x\|>t}\left|f^{*}(x)-f^{*}_{\Theta}(x/\|x\|)\right|\rightarrow
0\text{ as }t\rightarrow+\infty.$

The random pair $(X,Y)$ can be shown to fulfill Assumption 1 and conditional
regular variation of Assumption 2, refer to the Supplementary Material for
further details.

Here, by regression in the extremes is meant the usual regression problem in
regions far from the origin. More precisely, the objective pursued is the
construction of a real-valued, bounded and Borel-measurable function $f(x)$
minimizing the asymptotic conditional quadratic risk, defined as

$R_{\infty}(f):=\limsup_{t\to+\infty}R_{t}(f)=\limsup_{t\to+\infty}{\mathbb{E}}\left[(Y-f(X))^{2}\
|\ \|X\|\geq t\right].$ (12)

It is immediate to see that any function that coincides with
$f^{*}(x)={\mathbb{E}}\left[Y\mid X=x\right]$ on the region $\\{x\in
E,\|x\|>t\\}=t(E\setminus\mathbb{B})$ minimizes the risk functional $R_{t}$,
i.e. the quadratic risk taken w.r.t. distribution $P_{t}$. Set
$R_{t}^{*}:=R_{t}(f^{*})=\inf_{f}R_{t}(f)$ for all $t>0$ and observe that
$f^{*}$ also minimizes the asymptotic conditional quadratic risk functional:
$R_{\infty}^{*}:=R_{\infty}(f^{*})=\inf_{f}R_{\infty}(f)$.

As will be shown in the next section, the objective formulated above can be
connected to a classical regression problem related to the pair of random
variables $(X_{\infty},Y_{\infty})$ with distribution $P_{\infty}$ defined in
(10) by means of Assumptions 1 and 2 and describing the limit behavior of the
conditional distribution of the observable pair $(X/t,Y)$ given $\|X\|>t$ as
$t\to+\infty$. Precisely, given a class $\mathcal{F}$ of real-valued Borel
measurable and bounded functions defined on $E\setminus\mathbb{B}$, the goal
is to find $f$ in $\mathcal{F}$ that minimizes the quantity

$R_{P_{\infty}}(f):=\mathbb{E}\left[\left(Y_{\infty}-f(X_{\infty})\right)^{2}\right],$
(13)

referred to as the extreme quadratic risk. The best predictive performance one
can hope is achieved by the function $f_{P_{\infty}}^{*}$ defined by
$Y_{\infty}$’s conditional expectation given $X_{\infty}$,
$f_{P_{\infty}}^{*}(X_{\infty}):=\mathbb{E}[Y_{\infty}\mid X_{\infty}]$
almost-surely, and is denoted by
$R_{P_{\infty}}^{*}:=R_{P_{\infty}}(f_{P_{\infty}}^{*})=\inf_{f}R_{P_{\infty}}(f)$.
Based on its connection to extreme quadratic risk minimization, an algorithmic
approach to regression in the extremes can be designed. Before detailing it
and developing a sound theoretical framework to guarantee its validity, a few
remarks are in order.

###### Remark 2

(Heavy-tailed input vs heavy-tailed output) Attention should be paid to the
fact that the heavy-tail assumption is here on the distribution of the
input/explanatory r.v. $X$, in contrast to other works devoted to regression
such as Brownlees et al. (2015), Lugosi and Mendelson (2016) or Mendelson
(2017) where it is the loss/response that is supposedly heavy-tailed.

###### Remark 3

(Alternative to ERM) In the case where the output/response variable $Y$ is
heavy-tailed (or possibly contaminated by a heavy-tailed noise), alternatives
to the ERM approach Lugosi and Mendelson (2016). Extension of this approach is
beyond the scope of this paper but will be the subject of a forthcoming work.

##  3 Least Squares Regression in the Extremes

Under the assumptions previously listed, we now develop a rigorous framework
for least squares regression in extreme regions. We first establish
probabilistic results that connect the (asymptotic) quadratic conditional risk
and the extreme quadratic risk, as well as their minimizers. Based on the
latter, a statistical learning strategy to solve regression in the extremes is
next devised, a nonasymptotic upper confidence bound for the excess of
$R_{\infty}$-risk of the predictive function thus constructed.

###  3.1 Extreme Quadratic Risk Minimization through Regression in Extreme
Regions

Here, we shall establish the existence of angular predictive functions (i.e.
functions of the form $f(x)=(h\circ\theta)(x)$, where
$h:\mathbb{S}\to\mathbb{R}$ is a Borel measurable function) that
asymptotically minimize the $R_{\infty}$-risk, defined in (12). Based on this
result, we shall next propose a practical method to regression in the extremes
(see Algorithm 1) based on the angular information $\Theta(X_{i})$ provided by
a fraction of the largest observations $X_{i}$ (i.e. those with the largest
amplitude $||X_{i}||$) solely and prove its accuracy (see subsection 3.2).
Observe first that minimizers of the extreme quadratic risk (13) are of
angular type. Indeed, it follows from the definition (10) of the limit
distribution $P_{\infty}$ and the homogeneity of the $\mu_{y}$’s that the r.v.
$Y_{\infty}$ and $\|X_{\infty}\|$ are independent, so that:
$f_{P_{\infty}}^{*}(X_{\infty})={\mathbb{E}}\left[Y_{\infty}\mid
X_{\infty}\right]={\mathbb{E}}\left[Y_{\infty}\mid\Theta_{\infty}\right]$
almost-surely, where $\Theta_{\infty}:=\Theta(X_{\infty})$. Hence, the sole
part of information carried by $X_{\infty}$ that is useful to predict
$Y_{\infty}$ is its angular component $\Theta_{\infty}$. In order to express
the optimal solution $f_{P_{\infty}}^{*}$ in terms of the conditional angular
probability distributions, observe that $\Phi_{Y}$ is absolutely continuous
w.r.t. the marginal angular measure $\Phi$ with probability one and set
$\varphi_{Y}(\theta)=(d\Phi_{Y}/d\Phi)(\theta)$, for all
$\theta\in\mathbb{S}$.

###### Proposition 1

Suppose that Assumptions 1 and 2 are satisfied. For all $x\in E$, the
conditional distribution of $Y_{\infty}$ given $X_{\infty}=x$ (resp. given
$\Theta_{\infty}=\Theta(x)$) is absolutely continuous w.r.t. $G_{\infty}$ with
density

$g_{Y_{\infty}\mid
X_{\infty}=x}(y)=\frac{\varphi_{y}(\Theta(x))/\Phi_{y}(\mathbb{S})}{\int_{-M}^{M}\varphi_{u}(\Theta(x))/\Phi_{u}(\mathbb{S})dG_{\infty}(u)},\;\;y\in[-M,M].$

In addition, the minimizer of the extreme quadratic risk can be expressed as
$f_{P_{\infty}}^{*}(x)=f^{*}_{\Theta}(\Theta(x))$, where:
$\forall\theta\in\mathbb{S}$,

$f^{*}_{\Theta}(\theta)\ignorespaces=\frac{\int_{-M}^{M}u\varphi_{u}(\theta)/\Phi_{u}(\mathbb{S})dG_{\infty}(u)}{\int_{-M}^{M}\varphi_{u}(\theta)/\Phi_{u}(\mathbb{S})dG_{\infty}(u)}=\frac{\int_{-M}^{M}u\varphi_{u}(\theta)dG(u)}{\int_{-M}^{M}\varphi_{u}(\theta)dG(u)}.$

The proof mainly follows from the homogeneity property of $\mu_{y}$,
$y\in[-M,M]$, and is detailed in the Supplementary Material.

The conditional expectations and variances below are involved in the
following: $\forall x\in E$, $\forall\theta\in\mathbb{S}$,

$\displaystyle m_{P_{\infty}}(\theta)$
$\displaystyle={\mathbb{E}}\left[Y_{\infty}^{2}\mid\Theta_{\infty}=\theta\right],$
$\displaystyle V_{P_{\infty}}(\theta)$ $\displaystyle={\rm
Var}[Y_{\infty}\mid\Theta_{\infty}=\theta],$ $\displaystyle m(x)$
$\displaystyle={\mathbb{E}}\left[Y^{2}\mid X=x\right],$ $\displaystyle V(x)$
$\displaystyle={\rm Var}[Y\mid X=x].$

Equipped with these notations, the minimum extreme quadratic risk can be
expressed as follows:

$R_{P_{\infty}}^{*}=R_{P_{\infty}}(f_{P_{\infty}}^{*})=\mathbb{E}\Big{[}{\rm
Var}[Y_{\infty}\mid\Theta_{\infty}]\Big{]}.$

Observe also that: $\forall x\in E$, $\forall\theta\in\mathbb{S}$,

$V_{P_{\infty}}(\theta)=m_{P_{\infty}}(\theta)-(f^{*}_{\Theta}(\theta))^{2},\hskip
71.13188ptV(x)=m(x)-\left(f^{*}(x)\right)^{2}.$

The (mild) technical hypothesis below is also required.

###### Assumption 3

The functions $f^{*}_{\Theta}$ and $m_{P_{\infty}}$ are continuous on
$\mathbb{S}$. In addition, we have, as $t$ tends to infinity,

$\sup_{\theta\in\mathbb{S}}|f^{*}(t\theta)-f^{*}_{\Theta}(\theta)\mid\rightarrow
0\mbox{ and
}\sup_{\theta\in\mathbb{S}}|m(t\theta)-m_{P_{\infty}}(\theta)|\rightarrow 0.$
(14)

We point out that Assumption 3 is not very restrictive. Indeed, it is
automatically fulfilled as soon as additional classic mild conditions are
satisfied by the conditional density functions $\varphi_{y}$ De Haan and
Resnick (1987). Precisely, if for all $y\in[-M,M]$ the conditional
distribution of $X$ given $Y=y$ is absolutely continuous w.r.t. Lebesgue
measure on $\mathbb{R}^{d}$ with density $g_{X\mid Y=y}$ such that

$\sup_{\theta\in\mathbb{S}}|b(t)t^{d}g_{X\mid
Y=y}(t\theta)-\varphi_{y}(\theta)|\underset{t\rightarrow+\infty}{\longrightarrow}0,$
(15)

then the uniform convergences in (14) hold true. Refer to Appendix D and
Appendix E in the Supplementary Material for further details. Obviously, under
Assumption 3, $V_{P_{\infty}}$ is also continuous on $\mathbb{S}$ and the
uniform convergence below holds true:

$\sup_{x\in E:\|x\|\geq
t}|V_{P_{\infty}}(\Theta(x))-V(x)|\underset{t\rightarrow+\infty}{\longrightarrow}0.$
(16)

By relating asymptotic conditional quadratic risk minimization to extreme
quadratic risk minimization, the next proposition suggests to restrict the
search of minimizers of the $R_{\infty}$-risk to a class of angular predictive
functions.

###### Proposition 2

Suppose that Assumptions 1, 2 and 3 are satisfied. Then, the following
assertions hold true.

  * (i)

We have, as $t\rightarrow+\infty$:
$R_{t}^{*}\underset{t\rightarrow+\infty}{\longrightarrow}R_{P_{\infty}}^{*}$
and, in particular, $R_{\infty}^{*}=R_{P_{\infty}}^{*}$.

  * (ii)

The regression function $f_{P_{\infty}}^{*}=f^{*}_{\Theta}\circ\Theta$
minimizes the asymptotic conditional quadratic risk:

$R_{\infty}^{*}=R_{\infty}(f_{P_{\infty}}^{*}).$

The proposition stated above is proved in the Supplementary Material. It
reveals that the solution $f_{P_{\infty}}^{*}$ of the extreme risk
minimization problem, which is of angular type, is also a minimizer of the
asymptotic conditional quadratic risk (and that the minima coincide). Beyond
the assertions in Proposition 2, for any continuous function
$h:\mathbb{S}\rightarrow\mathbb{R}$, we also have:
$\lim_{t\rightarrow+\infty}R_{t}(h)=R_{\infty}(h)=R_{P_{\infty}}(h)$, see
Appendix A in the Supplementary Material. The next result, relying on the same
argument, shows that the conditional quadratic risk of the angular predictive
function $f_{P_{\infty}}^{*}$ converges to the minimum asymptotic conditional
quadratic risk.

###### Corollary 1

Suppose that Assumptions 1, 2 and 3 are satisfied. We have
$R_{t}(f_{\Theta}^{*}\circ\Theta)\underset{t\rightarrow+\infty}{\longrightarrow}R_{\infty}^{*}$
and, consequently,

$0\leq\inf_{h}R_{t}(h\circ\Theta)-R^{*}_{t}\rightarrow 0\text{ as
}t\rightarrow+\infty,$

where the infimum is taken over the class of all bounded Borel measurable
functions $h:\mathbb{S}\to\mathbb{R}$.

The corollary above encourages us to replace the original minimization problem
$\min R_{\infty}(f)$ by

$\min_{h}R_{t}(h\circ\Theta)$ (17)

with a threshold $t$ large enough, the bias resulting from the restriction to
the class of angular predictive functions being asymptotically negligible as
$t$ tends to $+\infty$. For $t$ suitably chosen (i.e. an appropriate order
statistic related to the $||X_{i}||$’s), an empirical version of the problem
(17) can be formulated. A practical and generic approach to regression in the
extremes then naturally consists in solving it, as described by Algorithm 1.

Algorithm 1 Least Squares Minimization in the Extremes

Input: Training data set
$\mathcal{D}_{n}=\\{(X_{1},Y_{1}),\;\ldots,\;(X_{n},Y_{n})\\}$ with
$X_{i}=(X_{i}^{(1)},\;\ldots,\;X_{i}^{(d)})\in[0,+\infty[^{d}$,
$i=1,\;\ldots,\;n$; class $\mathcal{H}$ of predictive functions on
$\mathbb{S}$; number $k\leq n$ of ’extreme’ observations selected.

Standardization: Standardize the input vectors by means of the transformation
$\hat{V}_{i}=\hat{T}(X_{i})\;\text{ for }i=1,\;\ldots,\;n,$ where
$\hat{T}(x)=\big{(}1/(1-\hat{F}_{1}(x^{(1)})),\;\ldots,\;1/(1-\hat{F}_{d}(x^{(d)}))\big{)}$,
with $\hat{F}_{j}(x^{(j)})=(1/(n+1))\sum_{i=1}^{n}\mathds{1}\\{X_{i}^{(j)}\leq
x^{(j)}\\}$ for $1\leq j\leq d$ and all $x\in\mathbb{R}^{d}$.

Truncation: Sort the training data by decreasing order of magnitude of the
transformed input $||\hat{V}_{(1)}||\geq\ldots\geq||\hat{V}_{(n)}||$ and form
a set of $k$ extreme training observations
$\left\\{\left(\hat{V}_{(1)},Y_{(1)}\right),\;\ldots,\;\left(\hat{V}_{(k)},Y_{(k)}\right)\right\\}.$
(18)

Empirical quadratic risk minimization: Solve the optimization problem
$\min_{h\in\mathcal{H}}\frac{1}{k}\sum_{i=1}^{k}\left(Y_{(i)}-h\left(\Theta\left(\hat{V}_{(i)}\right)\right)\right)^{2},$
(19) producing the solution $\hat{h}(\theta)$.

Output: Predictive function $(\hat{h}\circ\Theta)(x)$.

Notice that any algorithm for quadratic risk minimization can be used to solve
(19), refer to e.g. Györfi et al. (2002).

###### Remark 4

(Handling undefined quantities) In Algorithm 1, the empirical estimator of the
cdf $F_{j}$ is not exactly the usual version, the factor $1/n$ is replaced by
$1/(n+1)$ in order to avoid dividing by zero possibly in the marginal
transformation.

###### Remark 5

(Choosing the fraction $k$ of extreme observations) Choosing the number of
observations of a sample considered as extreme, i.e. the hyper-parameter $k$
in Algorithm 1, is a recurrent issue in EVT and a variety of methods to select
the optimal value of $k$ have been proposed in the literature, see e.g. Goix
et al. (2016), Goix et al. (2017) or Nakamura (2021) among others. However,
rules of thumb are generally used in practice, with the default choice
$k=\lfloor\sqrt{n}\rfloor$ in particular.

In the next subsection, a nonasymptotic analysis of the performance of the
approach proposed for regression in the extremes is carried out. An upper
confidence bound for the excess of $R_{\infty}$-risk of a solution of (19) is
established, when the class $\mathcal{H}$ over which empirical minimization is
performed is of controlled complexity. As previously mentioned, because the
goal of this paper is to explain the main ideas to tackle the problem of
regression in the extremes, for simplicity the input $1$-d marginal
distributions $F_{j}$ are supposedly unit Pareto. Hence, in order to avoid a
tedious technical analysis, the effect of step 1 of Algorithm 1 (empirical
standardization) is thus neglected here. It is however possible to incorporate
it in the rate bound analysis by means of the concentration results in
Clémençon et al. (2021) at the price of a significant amount of additional
technicalities, see the dedicated discussion in the Supplementary Material.

###  3.2 Theoretical Guarantees - Generalization Bounds

The rationale behind the approach proposed above to find a predictive function
that nearly minimizes the asymptotic conditional quadratic risk $R_{\infty}$
(12) consists in solving an empirical version of the non asymptotic
optimization problem

$\min_{h\in\mathcal{H}}R_{t}\left(h\circ\Theta\right),$

where the minimum is taken over a class $\mathcal{H}$ of continuous bounded
Borel measurable functions on $\mathbb{S}$ of controlled complexity (see
Assumption 4 below) but hopefully rich enough to contain a reasonable
approximant of $f^{*}_{\Theta}$ and $t$ is a large threshold. Large being
naturally meant in comparison with the observations of the training dataset
$\mathcal{D}_{n}$, the latter are first sorted by decreasing order of
magnitude of the $||X_{i}||$’s, $||X_{(1)}||>\ldots>||X_{(n)}||$, and one next
sets $t=\hat{t}_{n,k}:=||X_{(k)}||$. The threshold involved in the empirical
risk minimization problem is the statistical counterpart of $t_{n,k}$, the
quantile of level $1-k/n$ of $\|X\|$’s distribution. The empirical version of
the $R_{t_{n,k}}$-risk of a predictive mapping of the form $h\circ\Theta$ is

$\hat{R}_{\hat{t}_{n,k}}(h\circ\Theta):=\frac{1}{k}\sum_{i=1}^{k}\Big{(}Y_{(i)}-h(\Theta(X_{(i)}))\Big{)}^{2}.$
(20)

We point out that the statistic above is not an average of independent random
variables and that investigating its concentration properties is far from
straightforward.

###### Assumption 4

The class $\mathcal{H}$ is a set of continuous real-valued functions on
$\mathbb{S}$ of VC dimension $V_{\mathcal{H}}<+\infty$, uniformly bounded:
$\exists M<+\infty$ s.t. $\forall h\in\mathcal{H}$,
$\forall\theta\in\mathbb{S}$, $|h(\theta)|\leq M$.

Note that ’VC dimension’ refers here to the generalization of the usual VC
dimension for collections of sets to classes of real-valued functions, also
called Pollard’s pseudo-dimension sometimes Pollard (1990). Under the
complexity hypothesis above, the following result provides an upper confidence
bound for the maximal deviations between the conditional quadratic risk
$R_{t_{n,k}}$ and its empirical version, uniformly over the class
$\mathcal{H}$. It also shows that the bias inherent in the approximation of
the limit risk by a nonasymptotic version $R_{t}$ vanishes as
$t\rightarrow+\infty$, provided that the density of the conditional
distribution of $\Theta(X)$ given $\|X\|\geq t$, denoted by $\phi_{t}$, w.r.t.
$\Phi_{\Theta(X)}$, the marginal distribution of $\Theta(X)$, is uniformly
bounded. Refer to the Supplementary Material for its technical proof.

###### Proposition 3

Suppose that Assumptions 1-4 hold true.

  * (i)

Let $\delta\in(0,1)$. We have with probability larger than $1-\delta$:

$\sup_{h\in\mathcal{H}}\left|\hat{R}_{\hat{t}_{n,k}}(h\circ\Theta)-R_{t_{n,k}}(h\circ\Theta)\right|\leq\frac{4M^{2}}{\sqrt{k}}\Big{(}C\sqrt{V_{\mathcal{H}}}+2\sqrt{2\log(3/\delta)}\Big{)}+\frac{8M^{2}\log(3/\delta)}{3k},$

where $C$ is a universal constant.

  * (ii)

Suppose that $\sup_{t\geq
1,\;\theta\in\mathbb{S}}\phi_{t}(\theta):=D<+\infty$. We have:

$\sup_{h\in\mathcal{H}}\left|R_{t}(h\circ\Theta)-R_{\infty}(h\circ\Theta)\right|\rightarrow
0\text{ as }t\rightarrow+\infty.$

The corollary below provides an upper confidence bound for the excess of
$R_{\infty}$-risk of any solution $\hat{f}_{\Theta,k}$ of the problem

$\min_{h\in\mathcal{H}}\hat{R}_{\hat{t}_{n,k}}(h\circ\Theta).$

It immediately follows from the obvious bound

$\displaystyle R_{\infty}(\hat{f}_{\Theta,k}\circ\Theta)-R^{*}_{\infty}\leq$
$\displaystyle
2\sup_{h\in\mathcal{H}}|\hat{R}_{\hat{t}_{n,k}}(h\circ\Theta)-R_{t_{n,k}}(h\circ\Theta)|+\left(\inf_{h\in\mathcal{H}}R_{t_{n,k}}(h\circ\Theta)-R_{t_{n,k}}^{*}\right)+\left(R_{t_{n,k}}^{*}-R^{*}_{\infty}\right)$
$\displaystyle+R_{\infty}(\hat{f}_{\Theta,k}\circ\Theta)-R_{t_{n,k}}(\hat{f}_{\Theta,k}\circ\Theta),$

combined with Proposition 3.

###### Corollary 2

Suppose that Assumptions 1-4 are staisfied. Let $\delta\in(0,1)$. We have with
probability larger than $1-\delta$:

$\displaystyle R_{\infty}(\hat{f}_{\Theta,k}\circ\Theta)-R^{*}_{\infty}\leq$
$\displaystyle\frac{8M^{2}}{\sqrt{k}}\Big{(}C\sqrt{V_{\mathcal{H}}}+2\sqrt{2\log(3/\delta)}\Big{)}+\frac{16M^{2}\log(3/\delta)}{3k}+\left(\inf_{h\in\mathcal{H}}R_{t_{n,k}}(h\circ\Theta)-R_{t_{n,k}}^{*}\right)+\left(R_{t_{n,k}}^{*}-R^{*}_{\infty}\right)$
(21)
$\displaystyle+\sup_{h\in\mathcal{H}}\left|R_{t_{n,k}}(h\circ\Theta)-R_{\infty}(h\circ\Theta)\right|.$
(22)

As expected, $k$ being the number of extreme observations used to implement
the learning procedure, the stochastic error term in the bound stated above is
of order $O(1/\sqrt{k})$. Two types of bias error are also involved: one is
due to the restriction to the class $\mathcal{H}$ (model bias), while the
other results from the substitution of the conditional quadratic risk for its
asymptotic limit. In particular, under the assumptions stipulated, as $n$ and
$k$ tend to $+\infty$ so that $k/n\rightarrow 0$, we have $t_{n,k}\sim
n/k\rightarrow+\infty$ and the last two terms on the right hand side of (21)
vanish. A quantification of their decay rate would require to extend the
conditional multivariate regular variation property introduced in subsection
2.2 and specify second-order conditions, as in de Haan and Resnick (1996).

##  4 Numerical Experiments

We now investigate the performance of the approach previously described and
theoretically analyzed for regression in the extremes from an empirical
perspective on several simulated and real datasets. The MSE in extreme regions
of angular regression functions output by specific implementations of
Algorithm 1 are compared to those of the classic regression functions, learned
in a standard fashion. As a first go, we present the simulation results for an
additive noise model with heavy-tailed design (see Example 1):
$Y=f^{*}(X)+\varepsilon_{0}$, where
$f^{*}(x)=\beta^{T}\Theta(x)(1-1/2\sqrt{\|x\|})$ for $x\in\mathbb{R}^{d}$ and
$\|\cdot\|=\|\cdot\|_{2}$, the design $X$ is generated according to a
multivariate logistic distribution with dependence parameter equal to $\xi=1$
(which means that extreme observations occur along the axes Stephenson
(2003)), the input $1$-d marginals are standard Pareto with shape parameter
$\alpha=3$ and $\varepsilon_{0}=Z_{0}\mathds{1}\\{|Z_{0}|\leq 1\\}$ is the
noise, where $Z_{0}$ is a centered Gaussian r.v. with standard deviation
$\sigma=0.1$. The simulated data are of dimension $d=10$, the size of the
training dataset is $n_{train}=10.000$, truncated so as to keep only the
$k_{train}=223$ ($=\lfloor\sqrt{n_{train}}\rfloor$) largest observations to
learn regression functions in the extremes. The size of the test dataset is
$n_{test}=100.000$ and the $k_{test}=316$ ($=\lfloor\sqrt{n_{test}}\rfloor$)
largest instances are used to evaluate predictive performance in the extremes.

Table 1: Average MSE (and standard deviation) for regression functions trained
using all observations, extreme observations and angles of extreme
observations, over $20$ independent replications of the dataset generated by
the additive/multiplicative noise models.

Datasets/Models |  Train on $X$ |  Train on $X\mid\|X\|$ large |  Train on $\Theta\mid\|X\|$ large
---|---|---|---
Add. model: OLS | $0.359{\scriptstyle\pm 0.182}$ | $0.0518{\scriptstyle\pm 0.0285}$ | $\mathbf{0.0034{\scriptstyle\pm 0.0003}}$
SVR | $0.618{\scriptstyle\pm 0.279}$ | $0.0627{\scriptstyle\pm 0.0429}$ | $\mathbf{0.0041{\scriptstyle\pm 0.0006}}$
tree | $0.059{\scriptstyle\pm 0.012}$ | $0.0267{\scriptstyle\pm 0.0066}$ | $\mathbf{0.0203{\scriptstyle\pm 0.0102}}$
RF | $0.053{\scriptstyle\pm 0.012}$ | $0.0189{\scriptstyle\pm 0.0051}$ | $\mathbf{0.0110{\scriptstyle\pm 0.0031}}$
Mult. model: OLS | $1.72{\scriptstyle\pm 4.40}$ | $0.1402{\scriptstyle\pm 0.2462}$ | $\mathbf{0.0029{\scriptstyle\pm 0.0005}}$
SVR | $0.0091{\scriptstyle\pm 0.0020}$ | $0.0072{\scriptstyle\pm 0.0017}$ | $\mathbf{0.0067{\scriptstyle\pm 0.0006}}$
tree | $0.014{\scriptstyle\pm 0.004}$ | $0.011{\scriptstyle\pm 0.004}$ | $\mathbf{0.003{\scriptstyle\pm 0.002}}$
RF | $0.008{\scriptstyle\pm 0.001}$ | $0.006{\scriptstyle\pm 0.001}$ | $\mathbf{0.002{\scriptstyle\pm 0.001}}$

We implemented four different regression algorithms by means of scikit-learn
(Pedregosa et al., 2011) with the default parameters, namely Ordinary Least
Squares (OLS), Support Vector Regression (SVR) (with linear kernel for
additive noise and Gaussian kernel for multiplicative noise), CART (tree) and
Random Forest (RF). Predictive functions have been learned using the full
training dataset, the truncated version composed of the $k_{train}$ largest
observations and the angles of the truncated version: $e=20$ independent
replications of the experiment have been performed. Table 1 displays the
performance of the algorithms considered for the simulated additive noise
model. A similar experiment has been carried out for the multiplicative noise
model with heavy-tailed design $Y=\varepsilon_{1}f^{*}(X)$, where
$f^{*}(x)=\mbox{cos}(1/\|x\|)\sum_{i=1}^{d/2}\Theta(x)_{2i-1}\mbox{sin}(\Theta(x)_{2i}\pi)$,
the design $X$ is generated according to a multivariate logistic distribution
with dependence parameter equal to $\xi=0.7$, the input 1-d marginals are
standard Pareto with shape parameter $\alpha=2$ and
$\varepsilon_{1}=Z_{1}\mathds{1}\\{0\leq Z_{1}\leq 2\\}$ is the noise, $Z_{1}$
being a Gaussian r.v. with mean $\mu=1$ and standard deviation $\sigma=0.1$.
We consider the case of dimension $d=16$, with $n_{train}=10.000$ and
$n_{test}=100.000$ and then, $k_{train}=223$ and $k_{test}=316$, and $e=20$.
The averaged MSE are displayed in Table 1. The latter shows that, in both
experiments, the approach we promote for regression in the extremes and
analyzed in the previous section clearly outperforms its competitors, no
matter the algorithm (i.e. the model bias) considered. This paper being the
first to consider regression in the extremes (see Remark 3 for a description
of regression problems of different nature with heavy-tailed data), no other
alternative approach is documented in the literature. Encouraged by this first
agreement between theoretical and numerical results, experiments on real data
have been conducted. Four datasets are considered: dataset_sales ($10738$
instances and $15$ attributes) from OpenML repository
(https://www.openml.org/), bank32NH ($8192$ instances and $32$ attributes)
from LIACC repository (https://www.dcc.fc.up.pt/ ltorgo/Regression), CCPP
(9568 instances and 4 attributes) and CASP (45730 instances and 9 attributes)
from UCI repository (https://archive.ics.uci.edu/). Each dataset has been
randomly split $20$ times into a test set of size $n_{test}$, with one third
of the original dataset, and a train set of size $n_{train}$ with the rest.
The number of training extreme observations is $k_{train}=\lfloor
0.1*n_{train}\rfloor$ and $k_{test}=\lfloor 0.1*n_{test}\rfloor$. Results are
summarized in Table 2. Our approach to regression in the extremes generally
surpasses standard regression or regression using the truncated sample, except
in a few situations: when the ’extreme regime’ is not fully reached,
regression based on $X\mid\|X\|$ large may outperform regression based on
$\Theta\mid\|X\|$ large (see the results for bank32NH) and the impact of the
model bias, reflected in the bound of Corollary 2 should not be ignored (see
the results for CASP, best for the SVR algorithm). Practical application of
the methodology proposed/analyzed here in high-dimensional situations is of
crucial importance, though beyond the scope of the present paper. As the
approach promoted can be combined with any regression method, it should be
then naturally coupled with appropriate algorithms, based on regularization
techniques in particular.

Table 2: Average MSE (and standard deviation) for predictive functions learned
using all observations, extremes and angles of the extreme observations, over
$20$ random splits of each of the 4 datasets.

Datasets/Models |  Train on $X$ |  Train on $X\mid\|X\|$ large |  Train on $\Theta\mid\|X\|$ large
---|---|---|---
dataset_sales: OLS | $65.9{\scriptstyle\pm 7.8}$ | $57.7{\scriptstyle\pm 6.0}$ | $\mathbf{39.1{\scriptstyle\pm 4.9}}$
SVR | $65.5{\scriptstyle\pm 10.0}$ | $70.5{\scriptstyle\pm 9.4}$ | $\mathbf{46.8{\scriptstyle\pm 6.2}}$
tree | $59.9{\scriptstyle\pm 8.7}$ | $\mathbf{37.8{\scriptstyle\pm 7.0}}$ | $46.1{\scriptstyle\pm 8.2}$
RF | $31.7{\scriptstyle\pm 3.1}$ | $\mathbf{24.2{\scriptstyle\pm 3.7}}$ | $27.4{\scriptstyle\pm 3.7}$
bank32NH: OLS | $0.019{\scriptstyle\pm 0.001}$ | $0.018{\scriptstyle\pm 0.002}$ | $\mathbf{0.014{\scriptstyle\pm 0.001}}$
SVR | $0.021{\scriptstyle\pm 0.001}$ | $0.015{\scriptstyle\pm 0.001}$ | $\mathbf{0.014{\scriptstyle\pm 0.001}}$
tree | $0.040{\scriptstyle\pm 0.003}$ | $\mathbf{0.021{\scriptstyle\pm 0.004}}$ | $0.022{\scriptstyle\pm 0.002}$
RF | $0.021{\scriptstyle\pm 0.001}$ | $\mathbf{0.010{\scriptstyle\pm 0.001}}$ | $0.011{\scriptstyle\pm 0.001}$
CCPP: OLS | $\mathbf{22.4{\scriptstyle\pm 3.1}}$ | $348.8{\scriptstyle\pm 13.9}$ | $59.5{\scriptstyle\pm 4.5}$
SVR | $310.8{\scriptstyle\pm 15.8}$ | $271.0{\scriptstyle\pm 14.8}$ | $\mathbf{60.9{\scriptstyle\pm 4.9}}$
tree | $28.3{\scriptstyle\pm 5.4}$ | $\mathbf{22.4{\scriptstyle\pm 4.5}}$ | $73.5{\scriptstyle\pm 12.3}$
RF | $15.1{\scriptstyle\pm 2.9}$ | $\mathbf{13.4{\scriptstyle\pm 2.7}}$ | $41.5{\scriptstyle\pm 4.4}$
CASP: OLS | $\mathbf{14.5{\scriptstyle\pm 1.3}}$ | $43.7{\scriptstyle\pm 1.2}$ | $15.7{\scriptstyle\pm 1.2}$
SVR | $13.5{\scriptstyle\pm 1.0}$ | $28.9{\scriptstyle\pm 1.3}$ | $\mathbf{11.9{\scriptstyle\pm 1.1}}$
tree | $\mathbf{14.2{\scriptstyle\pm 2.0}}$ | $14.5{\scriptstyle\pm 1.6}$ | $19.1{\scriptstyle\pm 2.0}$
RF | $8.0{\scriptstyle\pm 0.7}$ | $\mathbf{7.7{\scriptstyle\pm 0.8}}$ | $9.5{\scriptstyle\pm 1.1}$

##  5 Conclusion

We have provided a sound strategy to solve ’regression in the extremes’. The
asymptotic framework we have developed crucially relies on the (novel) notion
of conditional regular variation. When the distribution of the input $X$
conditioned upon the output $Y$ is regularly varying, the problem can be
stated and analyzed in a rigorous manner. We have described the optimal
solution and proved that it can be nearly recovered with nonasymptotic
guarantees by implementing a variant of the ERM principle, based on the
angular information carried by a fraction of the largest observations only. We
have also carried out numerical experiments to support the approach promoted,
highlighting the necessity of using a dedicated methodology to perform
regression in the extremes with guarantees.

## References

  * Györfi et al. [2002] László Györfi, Michael Kohler, Adam Krzyzak, and Harro Walk.  _A Distribution-Free Theory of Nonparametric Regression_.  Springer, 2002.
  * Massart [2007] Pascal Massart.  _Concentration Inequalities and Model Selection: Ecole d’Eté de Probabilités de Saint-Flour XXXIV, volume 1896 of Lecture Notes in Mathematics._ Springer-Verlag, 2007.
  * Lecué and Mendelson [2016] Guillaume Lecué and Shahar Mendelson.  _Learning subgaussian classes: Upper and minimax bounds_.  2016\.
  * Boucheron et al. [2013] Stéphane Boucheron, Gábor Lugosi, and Pascal Massart.  _Concentration Inequalities: A Nonasymptotic Theory of Independence_.  Oxford University Press, 2013.
  * Thomas et al. [2017] Albert Thomas, Stéphan Clemencon, Alexandre Gramfort, and Anne Sabourin.  Anomaly Detection in Extreme Regions via Empirical MV-sets on the Sphere.  In Aarti Singh and Jerry Zhu, editors, _Proceedings of the 20th International Conference on Artificial Intelligence and Statistics_ , volume 54 of _Proceedings of Machine Learning Research_ , pages 1011–1019, Fort Lauderdale, FL, USA, 20–22 Apr 2017. PMLR.  URL http://proceedings.mlr.press/v54/thomas17a.html.
  * Goix et al. [2017] Nicolas Goix, Anne Sabourin, and Stephan Clémençon.  Sparse representation of multivariate extremes with applications to anomaly detection.  _Journal of Multivariate Analysis_ , 161:12–31, 2017.
  * Jalalzai et al. [2018] Hamid Jalalzai, Stephan Clémençon, and Anne Sabourin.  On binary classification in extreme regions.  In _Advances in Neural Information Processing Systems_ , pages 3092–3100, 2018.
  * Resnick [2013] Sidney I. Resnick.  _Extreme values, regular variation and point processes_.  Springer, 2013.
  * Ohannessian and Dahleh [2012] Mesrob I. Ohannessian and Munther A. Dahleh.  Rare probability estimation under regularly varying heavy tails.  In _Conference on Learning Theory_ , pages 21–1, 2012.
  * Bubeck et al. [2013] Sébastien Bubeck, Nicolo Cesa-Bianchi, and Gábor Lugosi.  Bandits with heavy tail.  _IEEE Transactions on Information Theory_ , 59(11):7711–7717, 2013.
  * Carpentier and Valko [2014] Alexandra Carpentier and Michal Valko.  Extreme bandits.  In _Advances in Neural Information Processing Systems 27_ , pages 1089–1097. Curran Associates, Inc., 2014.
  * Achab et al. [2017] Mastane Achab, Stephan Clémençon, Aurélien Garivier, Anne Sabourin, and Claire Vernade.  Max k-armed bandit: On the ExtremeHunter algorithm and beyond.  In _Machine Learning and Knowledge Discovery in Databases_ , pages 389–404. Springer International Publishing, 2017.  doi:10.1007/978-3-319-71246-8_24.  URL https://doi.org/10.1007%2F978-3-319-71246-8_24.
  * Goix et al. [2015] Nicolas Goix, Anne Sabourin, and Stephan Clémençon.  Learning the dependence structure of rare events: a non-asymptotic study.  In _Conference on Learning Theory_ , pages 843–860, 2015.
  * Brownlees et al. [2015] Christian Brownlees, Emilien Joly, and Gábor Lugosi.  Empirical risk minimization for heavy-tailed losses.  _The Annals of Statistics_ , 43(6):2507–2536, 2015.
  * Clémençon et al. [2021] Stéphan Clémençon, Hamid Jalalzai, Anne Sabourin, and Johan Segers.  Concentration bounds for the empirical angular measure with statistical learning applications, 2021.
  * Lugosi and Mendelson [2016] Gabor Lugosi and Shahar Mendelson.  Risk minimization by median-of-means tournaments.  _arXiv preprint arXiv:1608.00757_ , 2016.
  * Mendelson [2017] Shahar Mendelson.  On aggregation for heavy-tailed classes.  _Probability Theory and Related Fields_ , 168(3-4):641–674, 2017.
  * De Haan and Resnick [1987] Laurens De Haan and Sidney I. Resnick.  On regular variation of probability densities.  _Stochastic processes and their applications_ , 25:83–93, 1987.
  * Goix et al. [2016] Nicolas Goix, Anne Sabourin, and Stephan Clémençon.  Sparse representation of multivariate extremes with applications to anomaly ranking.  In _Artificial Intelligence and Statistics_ , pages 75–83, 2016.
  * Nakamura [2021] Chikara Nakamura.  On choice of hyper-parameter in extreme value theory based on machine learning techniques.  _arXiv preprint arXiv:2107.06074_ , 2021.
  * Pollard [1990] David Pollard.  Empirical processes: theory and applications.  Ims, 1990.
  * de Haan and Resnick [1996] Laurens de Haan and Sidney I. Resnick.  Second-order regular variation and rates of convergence in extreme-value theory.  _Ann. Prob_ , pages 97–124, 1996.
  * Stephenson [2003] Alec Stephenson.  Simulating multivariate extreme value distributions of logistic type.  _Extremes_ , 6(1):49–59, 2003.
  * Pedregosa et al. [2011] Fabian Pedregosa, Gael Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, Jake Vanderplas, Alexandre Passos, David Cournapeau, Matthieu Brucher, Matthieu Perrot, and Edouard Duchesnay.  Scikit-learn: Machine learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830, 2011\.
  * Resnick [1987] Sidney I. Resnick.  _Extreme Values, Regular Variation and Point Processes_.  Springer Series in Operations Research and Financial Engineering, 1987\.
  * Aghbalou et al. [2022] Anass Aghbalou, Patrice Bertail, François Portier, and Anne Sabourin.  Cross validation for rare events.  _arXiv preprint arXiv:2202.00488_ , 2022.
  * van der Vaart and Wellner [1996] A. W. van der Vaart and Jon Wellner.  _Weak Convergence and Empirical Processes: With Applications to Statistics_.  Springer Series in Statistics. Springer-Verlag, New York, 1996\.  ISBN 978-0-387-94640-5.  doi:10.1007/978-1-4757-2545-2.
  * Cai et al. [2011] Juan-Juan Cai, John H.J. Einmahl, and Laurens De Haan.  Estimation of extreme risk regions under multivariate regular variation.  _The Annals of Statistics_ , pages 1803–1826, 2011.
  * Girard and Stupfler [2015] Stéphane Girard and Gilles Stupfler.  Extreme geometric quantiles in a multivariate regular variation framework.  _Extremes_ , 18(4):629–663, 2015.
  * Drees et al. [2000] Holger Drees, Sidney Resnick, and Laurens de Haan.  How to make a hill plot.  _The Annals of Statistics_ , 28(1):254–274, 2000\.

The Supplementary Material is structured as follows. In Appendix A, the
conditional regular variation property introduced in Assumption 2 is
formulated in alternative and equivalent ways. Appendix B establishes the
properties of the additive noise model example, as well as those of two
additional examples. The technical proofs of the results stated in Section 3
are provided in Appendix C. Finally, additional comments and experiments about
Assumption 3 and Equation (15) are respectively discussed in Appendix D and
Appendix E.

##  Appendix A Conditional Multivariate Regular Variation

Here, we give alternative formulations of the conditional regular variation
condition introduced in Assumption 2, following in the footsteps of the
classical multivariate regular variation conditions, see Section 5 in Resnick
[1987]. The proof is omitted insofar as they can be obtained by following line
by line the arguments in Resnick [1987], when applied to the conditional
distribution of $X$ given $Y=y$.

###### Theorem A.1

The following assertions are equivalent.

  * (i)

There exist a transition kernel from $[-M,M]$ to $E$ given by
$(y,B)\in[-M,M]\times\mathcal{B}(E)\mapsto\mu_{y}(B)$, with $\mu_{y}$ a non-
null Radon measure on $\mathcal{B}(E)$ for every $y\in[-M,M]$, and a regularly
varying function $b$ with tail index $1$ such that

$\lim_{t\rightarrow+\infty}b(t)\mathbb{P}\left\\{t^{-1}X\in B\mid
Y=y\right\\}=\mu_{y}(B)$

for all $y\in[-M,M]$, for all $B\in\mathcal{B}(E)$ bounded away from zero and
s.t. $\mu_{y}(\partial B)=0$.

  * (ii)

There exist a transition kernel from $[-M,M]$ to $E$ given by
$(y,B)\in[-M,M]\times\mathcal{B}(E)\mapsto\mu_{y}(B)$, with $\mu_{y}$ a non-
null Radon measure on $\mathcal{B}(E)$ for every $y\in[-M,M]$, and a regularly
varying function $b$ with tail index $1$ for any bounded and continuous
function $h:E\rightarrow\mathbb{R}$ that vanishes in a neighborhood of zero

$\lim_{t\rightarrow+\infty}b(t)\mathbb{E}\left[h(t^{-1}X)\mid
Y=y\right]=\int_{x\in E}h(x)d\mu_{y}(x)$

for all $y\in[-M,M]$.

  * (iii)

There exists a transition kernel from $[-M,M]$ to $\mathbb{S}$ given by
$(y,A)\in[-M,M]\times\mathcal{B}(\mathbb{S})\mapsto\Phi_{y}(A)$, with
$\Phi_{y}$ a non-null finite measure on $\mathcal{B}(\mathbb{S})$ for all
$y\in[-M,M]$, such that for all $r\geq 1$ and $A\in\mathcal{B}(\mathbb{S})$,

$\frac{\mathbb{P}\\{\Theta(X)\in A,\|X\|\geq tr\mid
Y=y\\}}{\mathbb{P}\\{\|X\|\geq t\mid
Y=y\\}}\underset{t\rightarrow+\infty}{\longrightarrow}c_{y}r^{-1}\Phi_{y}(A)$

for all $y\in[-M,M]$, where $c_{y}=\Phi_{y}(\mathbb{S})^{-1}$.

  * (iv)

There exists a transition kernel from $[-M,M]$ to $\mathbb{S}$ given by
$(y,A)\in[-M,M]\times\mathcal{B}(\mathbb{S})\mapsto\Phi_{y}(A)$, with
$\Phi_{y}$ a non-null finite measure on $\mathcal{B}(\mathbb{S})$ for all
$y\in[-M,M]$, such that for all $A\in\mathcal{B}(\mathbb{S})$,

$\mathbb{P}\\{\Theta(X)\in A\mid\|X\|\geq
t,Y=y\\}\underset{t\rightarrow+\infty}{\longrightarrow}c_{y}\Phi_{y}(A)$

for all $y\in[-M,M]$, where $c_{y}=\Phi_{y}(\mathbb{S})^{-1}$ and also for all
$r\geq 1$,

$\frac{\mathbb{P}\\{\|X\|\geq tr\mid Y=y\\}}{\mathbb{P}\\{\|X\|\geq t\mid
Y=y\\}}\underset{t\rightarrow+\infty}{\longrightarrow}r^{-1}$

for all $y\in[-M,M]$.

In particular, with those notations,

$\mu_{y}(\\{x\in E,\|x\|\geq r,\Theta(x)\in A\\})=r^{-1}\Phi_{y}(A)$

for all $y\in[-M,M],r\geq 1$ and $A\in\mathcal{B}(\mathbb{S})$.

##  Appendix B Examples

In this section, we show that the additive noise model (Example 1) satisfies
the requirements of Assumption 2. In addition, we develop an extra example,
namely Multiplicative Noise Model.

###  B.1 Additive Noise Model with Heavy-Tailed Random Design

The next proposition shows that the typical example of an additive noise model
is compatible with our framework. We provide explicit expressions for the
limit measures and multiplicative constants.

###### Proposition B.1

Consider the additive noise model (11)

$Y=f^{*}(X)+\varepsilon_{0},$

where $X$ is a random vector with distribution $\rho$ taking its values in
$[0,\infty)^{d}$ satisfying Property (4) with associated exponent measure
$\mu$ and scaling function $b$, $\varepsilon_{0}$ is a bounded and centered
r.v. defined on the same probability space independent from $X$, drawn from a
continuous distribution with continuous density function $g_{0}$ w.r.t. the
Lebesgue measure and $f^{*}:[0,\infty)^{d}\rightarrow\mathbb{R}$ is a Borelian
bounded function s.t. there exists a continuous mapping
$f_{\Theta}^{*}:\mathbb{S}\rightarrow\mathbb{R}$ fulfilling the condition

$\sup_{||x||>t}\left|f^{*}(x)-f_{\Theta}^{*}(x/||x||)\right|\rightarrow
0\text{ as }t\rightarrow+\infty.$ (23)

Then, the random pair $(X,Y)$ satisfies Assumption 1 with
$M=\sup_{x\in[0,\infty)^{d}}|f^{*}(x)|+||\varepsilon_{0}||_{\infty}$ and
conditional regular variation of Assumption 2 with transition kernel defined
by:

$\forall
y\in[-M,M],\;\;d\mu_{y}(x)=\kappa_{y}^{-1}g_{0}(y-f_{\Theta}^{*}(\Theta(x)))d\mu(x),$

where $\kappa_{y}=\int_{z\in E}g_{0}(y-f^{*}(z))d\rho(z)$.

proof. First, Assumption 1 is obviously fulfilled with
$M=\sup_{x\in[0,\infty)^{d}}|f^{*}(x)|+||\varepsilon_{0}||_{\infty}$, with
$\|\varepsilon_{0}\|_{\infty}$ the essential supremum of the bounded r.v.
$\varepsilon_{0}$. Second, for any bounded and continuous function
$l:E\rightarrow\mathbb{R}$ that vanishes in a neighborhood $V$ of the origin,
we have for all $y\in[-M,M]$ and $t>0$,

$b(t)\mathbb{E}\left[l(t^{-1}X)\mid Y=y\right]=b(t)\int_{x\in
E}l(t^{-1}x)g_{0}(y-f^{*}(x))d\rho(x)/\kappa_{y}.$

We have the decomposition:

$b(t)\int_{x\in E}l(t^{-1}x)g_{0}(y-f^{*}(x))d\rho(x)=b(t)\int_{x\in
E}l(t^{-1}x)g_{0}(y-f_{\Theta}^{*}(\Theta(x)))d\rho(x)\\\ +b(t)\int_{x\in
t(E\setminus
V)}l(t^{-1}x)\Big{(}g_{0}(y-f^{*}(x))-g_{0}(y-f_{\Theta}^{*}(\Theta(x)))\Big{)}d\rho(x).$
(24)

Observe that the second term on the right hand side of the equation above is
bounded in amplitude by

$b(t)\int_{x\in E}|l(t^{-1}x)|d\rho(x)\sup_{x\in t(E\setminus
V)}\left|g_{0}(y-f^{*}(x))-g_{0}(y-f_{\Theta}^{*}(\Theta(x)))\right|.$

As $l$ is a bounded and continuous function that vanishes in the neighborhood
$V$ of zero, we have by virtue of the regular variation property of $X$’s
distribution (see (ii) of Theorem A.1)

$b(t)\int_{x\in E}|l(t^{-1}x)|d\rho(x)\rightarrow\int_{x\in E}|l(x)|d\mu(x)$

as $t\rightarrow+\infty$. As the density $g_{0}$ is (uniformly) continuous on
its compact support and

$\sup_{x\in t(E\setminus
V)}\left|f^{*}(x)-f_{\Theta}^{*}(\Theta(x))\right|\rightarrow 0$

as $t\rightarrow+\infty$, by virtue of Equation (23), we obtain that the
second term on the right hand side of (24) vanishes as $t$ tends to infinity.
In addition, define for all $y\in[-M,M]$ the mapping on $E$

$l_{y}(x)=l(x)g_{0}(y-f_{\Theta}^{*}(\Theta(x))).$

As, for all $y\in[-M,M]$, $l_{y}$ is a bounded and continuous function (since
$l,g_{0}$ and $f_{\Theta}^{*}$ are continuous and $g_{0}$ is compactly
supported) that vanishes in the neighborhood $V$ of zero, we have by virtue of
the regular variation property of $X$’s distribution:

$b(t)\int
l(t^{-1}x)g_{0}(y-f_{\Theta}^{*}(\Theta(x)))d\rho(x)\rightarrow\int_{x\in
E}l_{y}(x)d\mu(x),$

as $t\rightarrow+\infty$. The proof is finished by noticing that, for all
$y\in[-M,M]$, we have, again by (ii) of Theorem A.1:

$b(t)\mathbb{E}\left[l(t^{-1}X)\mid
Y=y\right]\underset{t\rightarrow+\infty}{\longrightarrow}\int_{x\in
E}l_{y}(t^{-1}x)d\mu(x)/\kappa_{y}=\int_{x\in E}l(x)d\mu_{y}(x).$

$\blacksquare$

###  B.2 Multiplicative Noise Model with Heavy-Tailed Random Design

The framework we consider here for regression in extreme regions also
encompasses multiplicative noise models, used in a wide range of applications,
such as radar imagery for instance. We have the following result.

###### Proposition B.2

Consider the multiplicative noise model

$Y=\varepsilon_{1}f^{*}(X),$

where $X$ is a random vector with distribution $\rho$ taking its values in
$[0,\infty)^{d}$ satisfying Property (4), $\varepsilon_{1}$ is a real-valued,
positive and bounded, r.v. defined on the same probability space, independent
from $X$, with mean $1$ and drawn from a continuous distribution with
continuous density $g_{1}$ w.r.t. the Lebesgue measure and
$f^{*}:[0,\infty)^{d}\rightarrow[1,+\infty[$ is a Borelian bounded function
s.t. there exists a continuous mapping
$f_{\Theta}^{*}:\mathbb{S}\rightarrow[1,+\infty[$ fulfilling the condition

$\sup_{||x||>t}\left|f^{*}(x)/f_{\Theta}^{*}(x/||x||)-1\right|\rightarrow
0\text{ as }t\rightarrow+\infty.$ (25)

Then, the random pair $(X,Y)$ satisfies Assumption 1 with
$M=\sup_{x\in[0,\infty)^{d}}|f^{*}(x)|\times||\varepsilon_{1}||_{\infty}$ and
conditional regular variation of Assumption 2 with transition kernel defined
by:

$\forall
y\in[-M,M],\;\;d\mu_{y}(x)=\xi_{y}^{-1}\frac{g_{1}(y/f_{\Theta}^{*}(\Theta(x)))}{f_{\Theta}^{*}(\Theta(x))}d\mu(x),$

where $\xi_{y}=\int_{z\in E}\frac{g_{1}(y/f^{*}(z))}{f^{*}(z)}d\rho(z)$.

proof. Observe first that Assumption 1 is satisfied with
$M=||\varepsilon_{1}||_{\infty}\sup_{x\in[0,\infty)^{d}}|f^{*}(x)|$, with
$\|\varepsilon_{1}\|_{\infty}$ the essential supremum of the bounded r.v.
$\varepsilon$. Next, notice that for any bounded and continuous function
$l:E\rightarrow\mathbb{R}$ that vanishes in a neighborhood $V$, we have for
all $y\in[-M,M]$ and $t>0$,

$b(t)\mathbb{E}\left[l(t^{-1}X)\mid Y=y\right]=b(t)\int_{x\in
E}l(t^{-1}x)\frac{g_{1}(y/f^{*}(x))}{f^{*}(x)}d\rho(x)/\xi_{y}.$

We have the decomposition:

$b(t)\int_{x\in
E}l(t^{-1}x)\frac{g_{1}(y/f^{*}(x))}{f^{*}(x)}d\rho(x)=b(t)\int_{x\in
E}l(t^{-1}x)\frac{g_{1}(y/f_{\Theta}^{*}(\Theta(x)))}{f_{\Theta}^{*}(\Theta(x))}d\rho(x)\\\
+b(t)\int_{x\in t(E\setminus
V)}l(t^{-1}x)\bigg{(}\frac{g_{1}(y/f^{*}(x))}{f^{*}(x)}-\frac{g_{1}(y/f_{\Theta}^{*}(\Theta(x)))}{f_{\Theta}^{*}(\Theta(x))}\bigg{)}d\rho(x).$
(26)

Define for all $y\in[-M,M]$ the mapping on $E$

$l_{y}(x)=l(x)\frac{g_{1}(y/f_{\Theta}^{*}(\Theta(x)))}{f_{\Theta}^{*}(\Theta(x))}.$

As, for all $y\in[-M,M]$, $l_{y}$ is a bounded and continuous function (since
$l,g_{1}$ and $f_{\Theta}^{*}$ are continuous and $g_{1}$ is compactly
supported) that vanishes in the neighborhood $V$ of zero, we have by virtue of
the regular variation property of $X$’s distribution (see (ii) of Theorem
A.1):

$b(t)\int_{x\in
E}\bigg{|}l(t^{-1}x)\frac{g_{1}(y/f_{\Theta}^{*}(\Theta(x)))}{f_{\Theta}^{*}(\Theta(x))}\bigg{|}d\rho(x)\rightarrow\int_{x\in
E}\bigg{|}l(x)\frac{g_{1}(y/f_{\Theta}^{*}(\Theta(x)))}{f_{\Theta}^{*}(\Theta(x))}\bigg{|}d\mu(x)=\int_{x\in
E}|l_{y}(x)|d\mu(x)$

as $t\rightarrow+\infty$.

Observe that the second term on the right hand side of the equation above is
bounded in amplitude by

$b(t)\int_{x\in
E}\bigg{|}l(t^{-1}x)\frac{g_{1}(y/f_{\Theta}^{*}(\Theta(x)))}{f_{\Theta}^{*}(\Theta(x))}\bigg{|}d\rho(x)\sup_{x\in
t(E\setminus
V)}\bigg{|}\frac{g_{1}(y/f^{*}(x))}{g_{1}(y/f_{\Theta}^{*}(\Theta(x)))}\frac{f_{\Theta}^{*}(\Theta(x))}{f^{*}(x)}-1\bigg{|}.$

As the density $g_{1}$ is (uniformly) continuous on its compact support and

$\sup_{x\in t(E\setminus
V)}\left|f^{*}(x)/f_{\Theta}^{*}(\Theta(x))-1\right|\rightarrow 0$

as $t\rightarrow+\infty$ by virtue of Equation (25), we obtain that the second
term on the right hand side of (26) vanishes as $t$ tends to infinity.

The proof is finished by noticing that, for all $y\in[-M,M]$, we have, again
by (ii) of Theorem A.1:

$b(t)\mathbb{E}\left[l(t^{-1}X)\mid
Y=y\right]\underset{t\rightarrow+\infty}{\longrightarrow}\int_{x\in
E}l_{y}(t^{-1}x)d\mu(x)/\xi_{y}=\int_{x\in E}l(x)d\mu_{y}(x).$

$\blacksquare$

###  B.3 Multiplicative and Additive Noise Model with Heavy-Tailed Random
Design

This last model is more general and covers the two other one. A proof of this
result immediately follows from the two other proofs and is thus omitted.

###### Proposition B.3

Consider the multiplicative and additive noise model

$Y=\varepsilon_{1}f^{*}(X)+\varepsilon_{0},$

where $X$ is a random vector with distribution $\rho$ taking its values in
$[0,\infty)^{d}$ satisfying Property (4), $\varepsilon_{1}$ is a real-valued,
positive and bounded, r.v. defined on the same probability space, independent
from $X$, with mean $1$ and drawn from a continuous distribution with
continuous density $g_{1}$ w.r.t. the Lebesgue measure, $\varepsilon_{0}$ is a
bounded and centered r.v. defined on the same probability space, independent
from $X$ and drawn from a continuous distribution with continuous density
$g_{0}$ w.r.t. the Lebesgue measure and
$f^{*}:[0,\infty)^{d}\rightarrow\mathbb{R}_{+}^{*}$ is a Borelian bounded
function s.t. there exists a continuous mapping
$f_{\Theta}^{*}:\mathbb{S}\rightarrow\mathbb{R}_{+}^{*}$ fulfilling the
condition

$\sup_{||x||>t}\left|f^{*}(x)-f_{\Theta}^{*}(x/||x||)\right|\rightarrow
0\text{ as }t\rightarrow+\infty.$

Then, the random pair $(X,Y)$ satisfies Assumption 1 with
$M=\sup_{x\in[0,\infty)^{d}}|f^{*}(x)|\times||\varepsilon_{1}||_{\infty}+||\varepsilon_{0}||_{\infty}$
and conditional regular variation of Assumption 2 with transition kernel
defined by:

$\forall
y\in[-M,M],\;\;d\mu_{y}(x)=\iota_{y}^{-1}\int_{z\in\mathbb{R}}g_{1}(z)g_{0}(y-zf_{\Theta}^{*}(\Theta(x)))dzd\mu(x)=\iota_{y}^{-1}\int_{z\in\mathbb{R}}\frac{g_{1}\Big{(}\frac{y-z}{f_{\Theta}^{*}(\Theta(x))}\Big{)}g_{0}(z)}{f_{\Theta}^{*}(\Theta(x))}dzd\mu(x),$

where

$\iota_{y}=\int_{x\in
E}\int_{z\in\mathbb{R}}g_{1}(z)g_{0}(y-zf^{*}(x))dzd\rho(x)=\int_{x\in
E}\int_{z\in\mathbb{R}}\frac{g_{1}\Big{(}\frac{y-z}{f^{*}(x)}\Big{)}g_{0}(z)}{f^{*}(x)}dzd\rho(x).$

##  Appendix C Technical Proofs

This section gathers technical proofs of Section 3. We recall for convenience
statements of results.

###  C.1 Proof of Proposition 1

###### Proposition 1

Suppose that Assumptions 1 and 2 are satisfied. For all $x\in E$, the
conditional distribution of $Y_{\infty}$ given $X_{\infty}=x$ (resp. given
$\Theta_{\infty}=\Theta(x)$) is absolutely continuous w.r.t. $G_{\infty}$ with
density

$g_{Y_{\infty}\mid
X_{\infty}=x}(y)=\frac{\varphi_{y}(\Theta(x))/\Phi_{y}(\mathbb{S})}{\int_{-M}^{M}\varphi_{u}(\Theta(x))/\Phi_{u}(\mathbb{S})dG_{\infty}(u)},\;\;y\in[-M,M].$

In addition, the minimizer of the extreme quadratic risk can be expressed as
$f_{P_{\infty}}^{*}(x)=f^{*}_{\Theta}(\Theta(x))$, where:
$\forall\theta\in\mathbb{S}$,

$f^{*}_{\Theta}(\theta)\ignorespaces=\frac{\int_{-M}^{M}u\varphi_{u}(\theta)/\Phi_{u}(\mathbb{S})dG_{\infty}(u)}{\int_{-M}^{M}\varphi_{u}(\theta)/\Phi_{u}(\mathbb{S})dG_{\infty}(u)}=\frac{\int_{-M}^{M}u\varphi_{u}(\theta)dG(u)}{\int_{-M}^{M}\varphi_{u}(\theta)dG(u)}.$

proof. The expression of $f_{P_{\infty}}^{*}$ immediately follows from the
existence and the expression of the conditional probability density
$g_{Y_{\infty}\mid X_{\infty}=x}$ of $Y_{\infty}$ given $X_{\infty}=x$ w.r.t.
$G_{\infty}$ since $f_{P_{\infty}}^{*}(x):={\mathbb{E}}\left[Y_{\infty}\mid
X_{\infty}=x\right]$. The proof is simply to show that $g_{Y_{\infty}\mid
X_{\infty}=x}$ exists for all $x\in E\backslash\mathbb{B}$ and is given by

$g_{Y_{\infty}\mid
X_{\infty}=x}(y)=\frac{\varphi_{y}(\Theta(x))/\Phi_{y}(\mathbb{S})}{\int_{-M}^{M}\varphi_{u}(\Theta(x))/\Phi_{u}(\mathbb{S})dG_{\infty}(u)},\forall
y\in[-M,M].$

Let $l_{1}:E\backslash\mathbb{B}\rightarrow\mathbb{R}$ be a bounded measurable
function and let $l_{2}:E\times[-M,M]\rightarrow\mathbb{R}$ be a non-negative
measurable function such that
${\mathbb{E}}\left[l_{2}(X_{\infty},Y_{\infty})\right]<+\infty$. By definition
of the conditional expectation, we want to prove that

${\mathbb{E}}\left[l_{1}(X_{\infty})l_{2}(X_{\infty},Y_{\infty})\right]={\mathbb{E}}\left[l_{1}(X_{\infty})\int_{-M}^{M}l_{2}(X_{\infty},u)g_{Y_{\infty}\mid
X_{\infty}}(u)dG_{\infty}(u)\right],$

which is sufficient to prove the existence of $g_{Y_{\infty}\mid
X_{\infty}=x}$ and the expression will straightforwardly follow. The existence
of a density $\varphi_{y}$ of $\Phi_{y}$ w.r.t. $\Phi$ implies that for all
$y\in[-M,M]$, $(R_{\infty},\Theta_{\infty})$ given $Y_{\infty}=y$ admits a
probability density (refer to Section 5.4.1. in Resnick [1987] for more
details) w.r.t. $Leb_{\geq 1}\otimes\Phi$ given by

$g_{(R_{\infty},\Theta_{\infty})\mid
Y_{\infty}=y}(r,\theta)=r^{-2}c_{y}\varphi_{y}(\theta),\forall y\in[-M,M].$

Thus, $(R_{\infty},\Theta_{\infty})$ also admits a probability density w.r.t.
$Leb_{\geq 1}\otimes\Phi$ given by

$g_{(R_{\infty},\Theta_{\infty})}(r,\theta)=r^{-2}\int_{-M}^{M}c_{u}\varphi_{u}(\theta)dG_{\infty}(u).$

Since ${\mathbb{E}}\left[l_{2}(X_{\infty},Y_{\infty})\right]<+\infty$, by
Fubini’s theorem

$\displaystyle{\mathbb{E}}\left[l_{1}(X_{\infty})l_{2}(X_{\infty},Y_{\infty})\right]$
$\displaystyle=\int_{y=-M}^{M}\int_{r=1}^{\infty}\int_{\theta\in\mathbb{S}}l_{1}(r\theta)l_{2}(r\theta,y)c_{y}\varphi_{y}(\theta)r^{-2}drd\Phi(\theta)dG_{\infty}(y)$
$\displaystyle=\int_{r=1}^{\infty}\int_{\theta\in\mathbb{S}}l_{1}(r\theta)\int_{y=-M}^{M}l_{2}(r\theta,y)c_{y}\varphi_{y}(\theta)dG_{\infty}(y)r^{-2}drd\Phi(\theta)$
$\displaystyle=\int_{r=1}^{\infty}\int_{\theta\in\mathbb{S}}l_{1}(r\theta)\int_{y=-M}^{M}l_{2}(r\theta,y)c_{y}\varphi_{y}(\theta)dG_{\infty}(y)\frac{\int_{-M}^{M}c_{u}\varphi_{u}(\theta)dG_{\infty}(u)}{\int_{-M}^{M}c_{u}\varphi_{u}(\theta)dG_{\infty}(u)}r^{-2}drd\Phi(\theta)$
$\displaystyle={\mathbb{E}}\left[l_{1}(X_{\infty})\int_{-M}^{M}l_{2}(X_{\infty},y)c_{y}\varphi_{y}(\Theta(X_{\infty}))\frac{dG_{\infty}(y)}{\int_{-M}^{M}c_{u}\varphi_{u}(\Theta(X_{\infty}))dG_{\infty}(u)}\right]$
$\displaystyle={\mathbb{E}}\left[l_{1}(X_{\infty})\int_{-M}^{M}l_{2}(X_{\infty},y)g_{Y_{\infty}\mid
X_{\infty}}(y)dG_{\infty}(y)\right],$

and the proof is down.  $\blacksquare$

###  C.2 Proof of Proposition 2

###### Proposition 2

Suppose that Assumptions 1, 2 and 3 are satisfied. Then, the following
assertions hold true.

  * (i)

We have, as $t\rightarrow+\infty$:
$R_{t}^{*}\underset{t\rightarrow+\infty}{\longrightarrow}R_{P_{\infty}}^{*}$
and, in particular, $R_{\infty}^{*}=R_{P_{\infty}}^{*}$.

  * (ii)

The regression function $f_{P_{\infty}}^{*}=f^{*}_{\Theta}\circ\Theta$
minimizes the asymptotic conditional quadratic risk:

$R_{\infty}^{*}=R_{\infty}(f_{P_{\infty}}^{*}).$

proof. (i) Since $R_{t}^{*}=R_{t}(f^{*})={\mathbb{E}}\left[V(X)\mid\|X\|\geq
t\right]$ and
$R_{P_{\infty}}^{*}={\mathbb{E}}\left[V_{P_{\infty}}(X_{\infty})\right]$, we
have

$\displaystyle|R_{t}^{*}-R_{P_{\infty}}^{*}|$
$\displaystyle=\Big{|}{\mathbb{E}}\left[V(X)\mid\|X\|\geq
t\right]-{\mathbb{E}}\left[V_{P_{\infty}}(X_{\infty})\right]\Big{|}$
$\displaystyle\leq\frac{\int_{\|x\|\geq
t}|V(x)-V_{P_{\infty}}(x)|d\rho(x)}{\mathbb{P}\\{\|X\|\geq
t\\}}+\bigg{|}\frac{\int_{\|x\|\geq
t}V_{P_{\infty}}(x)d\rho(x)}{\mathbb{P}\\{\|X\|\geq
t\\}}-{\mathbb{E}}\left[V_{P_{\infty}}(X_{\infty})\right]\bigg{|}$

The first term is controlled by Assumption 3:

$\frac{\int_{\|x\|\geq
t}|V(x)-V_{P_{\infty}}(x)|d\rho(x)}{\mathbb{P}\\{\|X\|\geq t\\}}\leq\sup_{x\in
E:\|x\|\geq t}|f^{*}(x)^{2}-f_{P_{\infty}}^{*}(x)^{2}|+\sup_{x\in E:\|x\|\geq
t}|m(x)-m_{P_{\infty}}(x)\mid\underset{t\rightarrow+\infty}{\longrightarrow}0.$

The second term tends to zero when $t$ tends to infinity since
$V_{P_{\infty}}$ is bounded and continuous and $X$ is regularly varying.
Finally, the equality $R_{\infty}^{*}=R_{P_{\infty}}^{*}$ ensues from the
beginning of the proof and
$R_{\infty}^{*}=\inf_{f}\limsup_{t}R_{t}(f)=\limsup_{t}R_{t}(f^{*})=\limsup
R_{t}^{*}$.

(ii) Notice that if a regression function $\tilde{f}$ satisfies

$\limsup_{t}\\{R_{t}(f^{*})-R_{t}(\tilde{f})\\}=0$ (27)

then $\tilde{f}$ minimizes $R_{\infty}$ as well. We shall prove that (27)
holds for $\tilde{f}=f_{P_{\infty}}^{*}$.

$\displaystyle R_{t}(f^{*})-R_{t}(f_{P_{\infty}}^{*})$
$\displaystyle={\mathbb{E}}\left[f_{P_{\infty}}^{*}(X)^{2}-f^{*}(X)^{2}\mid\|X\|\geq
t\right]+2{\mathbb{E}}\left[Y(f^{*}(X)-f_{P_{\infty}}^{*}(X))\mid\|X\|\geq
t\right]$ $\displaystyle\leq\sup_{\|x\|\geq
t}|f_{P_{\infty}}^{*}(x)^{2}-f^{*}(x)^{2}|+2M\sup_{\|x\|\geq
t}|f_{P_{\infty}}^{*}(x)-f^{*}(x)|$
$\displaystyle\underset{t\rightarrow+\infty}{\longrightarrow}0,$

when the last convergence follows from Assumption 3.  $\blacksquare$

###  C.3 Proof of Proposition 3

###### Proposition 3

Suppose that Assumptions 1-4 are satisfied.

  * (i)

Let $\delta\in(0,1)$. We have with probability larger than $1-\delta$:

$\sup_{h\in\mathcal{H}}\left|\hat{R}_{\hat{t}_{n,k}}(h\circ\Theta)-R_{t_{n,k}}(h\circ\Theta)\right|\leq\frac{4M^{2}}{\sqrt{k}}\Big{(}C\sqrt{V_{\mathcal{H}}}+2\sqrt{2\log(3/\delta)}\Big{)}+\frac{8M^{2}\log(3/\delta)}{3k},$

where $C$ is a universal constant.

  * (ii)

Let $\phi_{t}$ be the density of the conditional distribution of $\Theta(X)$
given $\|X\|\geq t$ w.r.t. $\Phi_{\Theta(X)}$, the marginal distribution of
$\Theta(X)$, and suppose that $\sup_{t\geq
1,\;\theta\in\mathbb{S}}\phi_{t}(\theta)=D<+\infty$. Then, as
$t\rightarrow+\infty$, we have:

$\sup_{h\in\mathcal{H}}\left|R_{t}(h\circ\Theta)-R_{\infty}(h\circ\Theta)\right|\rightarrow
0.$

proof. (i) The first assertion of Proposition 3, derives from Lemma A.3 of
Aghbalou et al. [2022], which we state below within our notational framework
for the sake of readability.

###### Lemma C.1

[Lemma A.3 of Aghbalou et al. [2022]] Let $\mathcal{H}$ be a family of
functions with finite VC-dimension $V_{\mathcal{H}}$. Then, the following
Bernstein inequality holds for all $t>0$

$\mathbb{P}\Big{\\{}\sup_{h\in\mathcal{H}}|\hat{R}_{t_{n,k}}(h\circ\Theta)-R_{t_{n,k}}(h\circ\Theta)|/(4M^{2})-\frac{C\sqrt{V_{\mathcal{H}}}}{\sqrt{k}}\geq
t\Big{\\}}\leq 3\exp\bigg{(}\frac{-kt^{2}}{2(4+t/3)}\bigg{)},$

where $C$ is an universal constant.

Assertion $(i)$ follows using the fact that $\forall a,b\geq
0,\sqrt{a+b}\leq\sqrt{a}+\sqrt{b}$.

(ii) We now prove the second assertion. Let $\phi_{\infty}$ be the density of
$\Theta_{\infty}$ w.r.t. $\Phi_{\Theta(X)}$. Observe first that for any
$h_{1},\;h_{2}$ in $\mathcal{H}$, we have:

$\left|R_{t}(h_{1}\circ\Theta)-R_{\infty}(h_{1}\circ\Theta)-\left(R_{t}(h_{2}\circ\Theta)-R_{\infty}(h_{2}\circ\Theta)\right)\right|\leq
4M\int_{\mathbb{S}}\left|h_{1}(\theta)-h_{2}(\theta)\right|\left(\phi_{t}(\theta)+\phi_{\infty}(\theta)\right)d\Phi_{\Theta(X)}(\theta)\\\
\leq
8MD\int_{\mathbb{S}}\left|h_{1}(\theta)-h_{2}(\theta)\right|d\Phi_{\Theta(X)}(\theta).$

Let $\varepsilon>0$ be fixed. Since $\mathcal{H}$ has a finite VC-dimension
(Assumption 4), the covering number of the class for the $L_{1}(\Phi)$-norm is
finite, $\mathcal{N}(\varepsilon,\;\mathcal{H},\;L_{1}(\Phi))<+\infty$ (see
e.g. van der Vaart and Wellner [1996], Section 2.6.2). Thus, for some
$L\leq\mathcal{N}(\varepsilon,\;\mathcal{H},\;L_{1}(\Phi))$, there exists
$h_{1},\;\ldots,\;h_{L}\in\mathcal{H}$ such that each $h\in\mathcal{H}$ is at
$L_{1}(\Phi)$-distance at most $\varepsilon$ from one of the $h_{i}$’s. In
other words for all $h\in\mathcal{H}$, $\exists i\leq L$ such that
$\|h-h_{i}\|_{L_{1}(\Phi)}\leq\varepsilon$. Thus by the triangular inequality

$\sup_{h\in\mathcal{H}}\left|R_{t}(h\circ\Theta)-R_{\infty}(h\circ\Theta)\right|\leq\max_{1\leq
l\leq
L}\left|R_{t}(h_{l}\circ\Theta)-R_{\infty}(h_{l}\circ\Theta)\right|+8MD\varepsilon.$

As the first term on the right hand side of the bound above is less than
$\varepsilon$ for $t$ large enough, the uniform convergence is proved. Notice
incidentally that the mild condition $\sup_{t\geq
1,\;\theta\in\mathbb{S}}\phi_{t}(\theta)=D<+\infty$ can be naturally weakened
and replaced by an appropriate moment condition.  $\blacksquare$

###  C.4 Proof of Corollary 2

###### Corollary 2

Suppose that Assumptions 1-4 are satisfied. Let $\delta\in(0,1)$. We have with
probability larger than $1-\delta$:

$\displaystyle R_{\infty}(\hat{f}_{\Theta,k}\circ\Theta)-R^{*}_{\infty}$
$\displaystyle\leq\frac{8M^{2}}{\sqrt{k}}\Big{(}C\sqrt{V_{\mathcal{H}}}+2\sqrt{2\log(3/\delta)}\Big{)}+\frac{16M^{2}\log(3/\delta)}{3k}+\left(\inf_{h\in\mathcal{H}}R_{t_{n,k}}(h\circ\Theta)-R_{t_{n,k}}^{*}\right)$
$\displaystyle+\left(R_{t_{n,k}}^{*}-R^{*}_{\infty}\right)+\sup_{h\in\mathcal{H}}\left|R_{t_{n,k}}(h\circ\Theta)-R_{\infty}(h\circ\Theta)\right|.$

proof. The proof relies on the inequality preceding Proposition 3 which itself
follows from the approximation/estimation risk decomposition given by

$|R_{t_{n,k}}(\hat{f}_{\Theta,k})-\inf_{f}R_{t_{n,k}}(f)|\leq
2\sup_{h\in\mathcal{H}}|\hat{R}_{t_{n,k}}(h\circ\Theta)-R_{t_{n,k}}(h\circ\Theta)|+\inf_{h\in\mathcal{H}}R_{t_{n,k}}(h\circ\Theta)-\inf_{f}R_{t_{n,k}}(f).$

Equipped with this inequality, Corollary 2 is just an application of
Proposition 3.  $\blacksquare$

##  Appendix D On Assumption 3 and Empirical Marginal Standardization

This uniform convergence assumption mentioned in the result below is stated in
De Haan and Resnick [1987] and used in Cai et al. [2011] or Girard and
Stupfler [2015]. A proof of the claim following Assumption 3 is detailed
below.

###### Proposition D.4

(Sufficient conditions for Assumption 3) Suppose that Assumptions 1 and 2 are
satisfied and suppose that, for all $y\in[-M,M],\Phi_{y}$ admits a probability
density function $\varphi_{y}$ w.r.t. $\lambda$, the Hausdorff surface measure
on $\mathbb{S}$. Suppose also that, for all $y\in[-M,M]$, there exists a
continuous probability density function of $X$ given $Y=y$ w.r.t. the Lebesgue
measure on $\mathbb{R}^{d}$, denoted by $g_{X\mid Y=y}$, such that for all
$y\in[-M,M]$,

$\sup_{\theta\in\mathbb{S}}|b(t)t^{d}g_{X\mid
Y=y}(t\theta)-\varphi_{y}(\theta)|\underset{t\rightarrow+\infty}{\longrightarrow}0.$

Then, for all $l:[-M,M]\rightarrow\mathbb{R}$ continuous, we have

$\sup_{\theta\in\mathbb{S}}\Big{|}{\mathbb{E}}\left[l(Y)\mid
X=t\theta\right]-{\mathbb{E}}\left[l(Y_{\infty})\mid\Theta_{\infty}=\theta\right]\Big{|}\underset{t\rightarrow+\infty}{\longrightarrow}0.$

proof. Let $l:[-M,M]\rightarrow\mathbb{R}$ be a continuous function. First,
similarly to maths performed in the proof of Proposition 1, we have

${\mathbb{E}}\left[l(Y)\mid X=t\theta\right]=\int_{-M}^{M}l(y)g_{X\mid
Y=y}(t\theta)\frac{dG(y)}{g_{X}(t\theta)}.$

Thanks to Proposition 1, one can write

${\mathbb{E}}\left[l(Y)\mid
X=t\theta\right]-{\mathbb{E}}\left[l(Y_{\infty})\mid\Theta_{\infty}=\theta\right]=\int_{-M}^{M}l(y)\bigg{(}\frac{g_{X\mid
Y=y}(t\theta)}{g_{X}(t\theta)}-\frac{\varphi_{y}(\theta)}{\int_{-M}^{M}\varphi_{u}(\theta)dG(u)}\bigg{)}dG(y).$

Observe that for all $y\in[-M,M]$

$\displaystyle\frac{g_{X\mid
Y=y}(t\theta)}{g_{X}(t\theta)}-\frac{\varphi_{y}(\theta)}{\int_{-M}^{M}\varphi_{u}(\theta)dG(u)}$
$\displaystyle=\frac{b(t)t^{d}g_{X\mid
Y=y}(t\theta)}{b(t)t^{d}g_{X}(t\theta)}-\frac{\varphi_{y}(\theta)}{\int_{-M}^{M}\varphi_{u}(\theta)dG(u)}$
$\displaystyle=\frac{1}{\int_{-M}^{M}\varphi_{u}(\theta)dG(u)}\Big{(}b(t)t^{d}g_{X\mid
Y=y}(t\theta)-\varphi_{y}(\theta)\Big{)}$ $\displaystyle+b(t)t^{d}g_{X\mid
Y=y}(t\theta)\bigg{(}\frac{1}{b(t)t^{d}g_{X}(t\theta)}-\frac{1}{\int_{-M}^{M}\varphi_{u}(\theta)dG(u)}\bigg{)}$

Then, it follows

$\sup_{\theta\in\mathbb{S}}\Big{|}{\mathbb{E}}\left[l(Y)\mid
X=t\theta\right]-{\mathbb{E}}\left[l(Y_{\infty})\mid\Theta_{\infty}=\theta\right]\Big{|}\leq\frac{1}{\min_{\theta\in\mathbb{S}}\Big{|}\int_{-M}^{M}\varphi_{u}(\theta)dG(u)\Big{|}}\int_{-M}^{M}|l(y)|\sup_{\theta\in\mathbb{S}}|b(t)t^{d}g_{X\mid
Y=y}(t\theta)-\varphi_{y}(\theta)|dG(y)\\\
+\frac{\int_{-M}^{M}\sup_{\theta\in\mathbb{S}}|b(t)t^{d}g_{X\mid
Y=u}(t\theta)-\varphi_{u}(\theta)|dG(u)}{\min_{\theta\in\mathbb{S}}|\
b(t)t^{d}g_{X}(t\theta)\int_{-M}^{M}\varphi_{u}(\theta)dG(u)|}\int_{-M}^{M}|l(y)|\sup_{\theta\in\mathbb{S}}|b(t)t^{d}g_{X\mid
Y=y}(t\theta)|dG(y).$

The existence of $\min_{\theta\in\mathbb{S}}|b(t)t^{d}g_{X}(t\theta)|$ follows
from the continuity of $g_{X\mid Y=y}$ and this quantity is different from
zero since
$\min_{\theta\in\mathbb{S}}\Big{|}\int_{-M}^{M}\varphi_{u}(\theta)dG(u)\Big{|}>0$
(a proof by contradiction merely shows this assertion). Since for all
$y\in[-M,M]$, $g_{X\mid Y=y}$ is continuous and Equation (15) holds,
$\sup_{t\geq 1}\sup_{\theta\in\mathbb{S}}|b(t)t^{d}g_{X\mid
Y=y}(t\theta)-\varphi_{y}(\theta)|$ and $\sup_{t\geq
1}\sup_{\theta\in\mathbb{S}}|b(t)t^{d}g_{X\mid Y=y}(t\theta)|$ are finite for
all $\theta\in\mathbb{S}$ and using a dominated convergence argument and the
uniform convergence of the densities (15), the proof is over.  $\blacksquare$

###### Remark 6

Under hypothesis of Proposition D.4, Assumption 3 is entirely fulfilled. The
continuity of $f_{P_{\infty}}^{*}$ and $m_{P_{\infty}}$ on $\mathbb{S}$
follows from continuity of $g_{X\mid Y=y}$ on $\mathbb{S}$ and a uniform
convergence argument (15). Then, taking successively $l:y\in[-M,M]\mapsto y$
and $l:y\in[-M,M]\mapsto y^{2}$ in Proposition D.4 involves that

$\sup_{\theta\in\mathbb{S}}|f^{*}(t\theta)-f_{\Theta}^{*}(\theta)\mid\rightarrow
0,\mbox{ and
}\sup_{\theta\in\mathbb{S}}|m(t\theta)-m_{P_{\infty}}(\theta)\mid\rightarrow
0,$

as $t\rightarrow+\infty$.

###### Remark 7

According to Theorem 1 of De Haan and Resnick [1987], condition (15) could be
a working hypothesis replacing our working hypothesis, i.e. Assumption 2.
Indeed, suppose that for all $y\in[-M,M]$ there exist a monotone continuous
probability density function $g_{X\mid Y=y}$ w.r.t. the Lebesgue measure on
$\mathbb{R}^{d}$, a limiting function $q_{y}$ and a regularly varying function
$b$ with tail index $1$ s.t. for all $x\in E$

$b(t)t^{d}g_{X\mid Y=y}(tx)\rightarrow q_{y}(x)$

as $t\rightarrow+\infty$. Then, the conditional distribution of $X$ given
$Y=y$ is regularly varying with exponent measure $\mu_{y}$ such that

$\mu_{y}(B)=\int_{B}q_{y}(x)dx.$

The reader is invited to consult De Haan and Resnick [1987] for further
details and examples.

###### Remark 8

(On empirical marginal standardization) In this paper, marginal distributions
are assumed to be known and, as stressed in Remark 1, technical issues arise
without this assumption. This is a recurrent theme in multivariate extreme
value theory. In concrete cases, an empirical marginal standardization is a
mandatory step in order to use our work, as described in Algorithm 1, while no
guarantees are given in the present paper due to space limitation. Proof
techniques are very similar to those used in the binary classification case,
where estimation error bounds remain of the order $O_{\mathbb{P}}(1/\sqrt{k})$
while working with empirically standardized observations (as proved in
Clémençon et al. [2021], Section 4.2), could be investigated. Such a
refinement is beyond the scope of this paper and will be the subject of
further research.

##  Appendix E Checking Moments Assumption: Empirical Methodology and Results

In this section we propose a simple heuristic aimed at checking our specific
regularity assumptions, namely the assumption of moments convergence,
Assumption 3. We investigate here the empirical ’pointwise’ convergence of the
empirical conditional regression function $f^{*}(t\cdot)$ above level $t$
towards the extreme regression function $f_{\Theta}^{*}$ on the sphere, with
the notations introduced in Section 3. A similar analysis may be carried out
regarding the convergence of second moments (second moment regression function
$m$), with similar results.

Since this work is intended to open the road to further studies regarding
regression problems in extreme regions we consider as a first go a simple
setting, namely a low-to-moderate dimensional case where histogram-based rules
are not suffering from the curse of dimensionality. The higher dimensional
case is left to future works. We thus consider a regular partition of the
sphere $\mathbb{S}$, denoted by $\\{\Omega_{1},...,\Omega_{p}\\}$ (e.g. the
intersection of a regular grid in $\mathbb{R}^{d}$ with the sphere, following
Clémençon et al. [2021]), and introduce smoothed versions of
$f^{*}(t\theta_{j})$ for some $\theta_{j}\in\Omega_{j}$,

$f_{t,j}:={\mathbb{E}}\left[Y\mid\|X\|\geq t,\Theta(X)\in\Omega_{j}\right],$

for every $j\in\\{1,...,p\\}$. As in Section 4, we consider in our experiments
the $L_{2}$-sphere in $\mathbb{R}^{d}$.

The proposed method consists in investigating the convergence of empirical
versions of the $f_{t,j}$’s as $t$ becomes large. Let
$((X_{1},Y_{1}),...,(X_{n},Y_{n}))$ be a sample of independent observations
drawn from the pair $(X,Y)\in\mathbb{R}^{d}\times\mathbb{R}$. Recall that the
order statistics sorted by decreasing order of magnitude are denoted by
$\|X_{(1)}\|>\ldots>\|X_{(n)}\|$ with $Y_{(i)}$ the corresponding sorted
labels. Notice that $\hat{t}_{n,k}:=\|X_{(k)}\|$ is the empirical version of
$t_{n,k}$, the quantile of level $1-k/n$ of $\|X\|$’s distribution.

We investigate the behavior of the empirical version of $f_{t_{n,k},j}$,

$\hat{f}_{\hat{t}_{n,k},j}:=\bigg{(}\sum_{i=1}^{n}Y_{i}\mathds{1}\\{\|X_{i}\|\geq\hat{t}_{n,k},\Theta_{i}\in\Omega_{j}\\}\bigg{)}\bigg{(}1+\sum_{i=1}^{n}\mathds{1}\\{\|X_{i}\|\geq\hat{t}_{n,k},\Theta_{i}\in\Omega_{j}\\}\bigg{)}^{-1},$

for $j\leq p$ and $k\leq n$. Notice that adding one to the sum of indicator
functions in the denominator is simply meant to prevent dividing by zero when
no observation belongs to a region $\Omega_{j}$. Under Assumption 3, for large
values of $n,k$ with $k\ll n$, the $\hat{f}_{\hat{t}_{n,k},j}$’s should
stabilize around the limit
$f_{\Theta,j}:=\mathbb{E}[Y_{\infty}\mid\Theta_{\infty}\in\Omega_{j}]$. From a
practical point of view, the graphs of the $\hat{f}_{\hat{t}_{n,k},j}$ as a
function of $k$ (for fixed, large $n$ and fixed $j$) should then exhibit a
stability region for moderate values of $k$, large enough for the law of large
numbers to apply but small enough so that $k/n$ is small, hence $t_{n,k}$ is a
large quantile.

We analyze the behavior of these estimators while $k$ grows but remains
sufficiently small compared to $n$, and the goal of these experiments is to
exhibit existence of stability regions corresponding to values of the
$f_{\Theta,j}$’s for $j\in\\{1,...,p\\}$. Observe that $k$ must be small
w.r.t. the size of the sample $n$ in order to ensure that the possible
stabilization occurs under an extreme regime; a stability region should appear
at $k=n$ due to the convergence of $\hat{f}_{\hat{t}_{n,k},j}$ to
$\mathbb{E}[Y\mid\Theta(X)\in\Omega_{j}]$ when $k$ tends to infinity but this
is of no interest for our purpose. The search for stability regions of
estimators based on the $k^{th}$ largest order statistics of a random sample
is a classical approach in extreme value analysis, in particular for extreme
value index estimation where they are commonly referred to as ‘Hill plots’,
see _e.g._ Drees et al. [2000] and the references therein. Once such stability
regions have been exhibited, they may serve as guidance in order to choose a
reasonable value of $k$ and achieve a satisfactory bias-variance trade-off.

The experiments reported in this section aim at verifying experimentally the
existence of such stability regions for $k$ as described above and at
investigating the potential improvement of the regression algorithms in terms
of mean squared error when $k$ is chosen in such a stability region. In our
first experiment we place ourselves in the setting of the additive noise model
as described in Section 4, in dimension $d=5$ and simulate a sample size
$n=20.000$. Figure 1 displays estimates $\hat{f}_{\hat{t}_{n,k},j}$ as
function of $k$, for those $j\in\\{1,...,p\\}$ such that a sufficient number
of extreme points are located within $\Omega_{j}$, namely, such that there
exists $k\in\\{1,...,n\\},\|X_{k}\|\geq\|X_{(10)}\|,\Theta_{k}\in\Omega_{j}$.
We denote by $\omega_{10}$ the set of indices $j$ verifying the latter
condition. In other words we consider only elements of the partition
containing at least one of the ten most extreme observations.

In a second experiment we consider a real dataset of moderate dimension used
in the experiments from Section 4, namely the CCPP dataset. The results are
gathered in Figure 2. In both experiments a stability region indeed exists,
roughly within the range $k\in[100,200]$ for simulated data (Figure 1) and
$k\in[200,1000]$ for the _CCPP_ dataset (Figure 2).

Finally we compare MSE-performances obtained in Section 4 with MSE-
performances obtained by selecting $k$, the number of observations considered
as extreme, in stability regions. We follow the experimental workflow
described in Section 4 with a different of $k_{train}=\lfloor
0.05*n_{train}\rfloor$ and $k_{test}=\lfloor 0.05*n_{test}\rfloor$. Table 3
displays MSE-performances for these sample sizes. Overall performances are
indeed improved compared with the results reported using the default values of
$k$, see Table 2 from the main paper.

[5pt]Plot for $k\leq n/5$.

[5pt]Plot for $k\leq n/20$.

Figure 1: Plots of $\hat{f}_{\hat{t}_{n,k},j}$ against $k$, where
$j\in\omega_{10}$ for a dataset generated from the additive noise model with
parameters $p=5,d=5,n=20.000,\beta=(0.04,0.96,0.94,0.45,0.13),\xi=1,\alpha=3$.

[5pt]Plot for $k\leq n/2$. [5pt]Plot for $k\leq n/10$.

Figure 2: Plots of $\hat{f}_{\hat{t}_{n,k},j}$ against $k$, where
$j\in\omega_{10}$ for the CCPP dataset. Table 3: Average MSE (and standard
deviation) for predictive functions learned using all observations, extreme
observations and angles of the extreme observations, over $20$ random splits
of CCPP dataset, with the number of extreme observations selected in stability
regions displayed in Figure 2.
Datasets |  Train on $X$ |  Train on $X\mid\|X\|$ large |  Train on $\Theta\mid\|X\|$ large
---|---|---|---
CCPP: OLS | $\mathbf{22.1{\scriptstyle\pm 3.4}}$ | $292.6{\scriptstyle\pm 20.2}$ | $51.1{\scriptstyle\pm 5.2}$
SVR | $331.7{\scriptstyle\pm 20.4}$ | $238.1{\scriptstyle\pm 23.1}$ | $\mathbf{61.5{\scriptstyle\pm 9.5}}$
tree | $25.0{\scriptstyle\pm 4.6}$ | $\mathbf{16.7{\scriptstyle\pm 3.4}}$ | $61.9{\scriptstyle\pm 12.3}$
RF | $13.3{\scriptstyle\pm 2.2}$ | $\mathbf{9.1{\scriptstyle\pm 1.8}}$ | $36.7{\scriptstyle\pm 6.6}$