# Spectral Learning of Multivariate Extremes

Marco Avella Medina  Department of Statistics
Columbia University<EMAIL_ADDRESS>,  Richard A. Davis  Department
of Statistics
Columbia University<EMAIL_ADDRESS>and  Gennady Samorodnitsky
School of Operations Research and Information Engineering
Cornell University<EMAIL_ADDRESS>
###### Abstract.

We propose a spectral clustering algorithm for analyzing the dependence
structure of multivariate extremes. More specifically, we focus on the
asymptotic dependence of multivariate extremes characterized by the angular or
spectral measure in extreme value theory. Our work studies the theoretical
performance of spectral clustering based on a random $k$-nearest neighbor
graph constructed from an extremal sample, i.e., the angular part of random
vectors for which the radius exceeds a large threshold. In particular, we
derive the asymptotic distribution of extremes arising from a linear factor
model and prove that, under certain conditions, spectral clustering can
consistently identify the clusters of extremes arising in this model.
Leveraging this result we propose a simple consistent estimation strategy for
learning the angular measure. Our theoretical findings are complemented with
numerical experiments illustrating the finite sample performance of our
methods.

###### Key words and phrases:

Angular measure, heavy tails, Laplacian, nearest neighbor graphs, regular
variation, spectral clustering

This research was partially supported by NSF grants DMS-2015379 (Avella Medina
and Davis) at Columbia and DMS-2015242 (Samorodnitsky) at Cornell.

##  1\. Introduction

Multivariate extremes arise when one or more of rare extreme events occur
simultaneously. They are of paramount importance for understanding
environmental risks such as fires or droughts since they are driven by joint
extremes of a number of meteorological variables. Similarly, catastrophic
financial events are also of a multivariate nature in financial systems driven
by core institutions that are connected. In the above examples one is
precisely interested in modeling the dependence between rare individual
extremes. Multivariate extreme value theory is an active research area that
provides tools for modeling such events.

The dependence structure between extreme observations can be complex and
typically characterized by different notions of dependence from the ones
arising in the non-extreme world. For this reason recent work has sought to
rethink various notions of sparsity for extremes (Goix et al., 2017; Meyer and
Wintenberger, 2019; Simpson et al., 2020), concentration inequalities (Goix et
al., 2015; Clémençon et al., 2021), conditional independence (Belkin and
Niyogi, 2003) and unsupervised learning (Chautru, 2015; Cooley and Thibaud,
2019; Janßen and Wan, 2020; Drees and Sabourin, 2021). See also Engelke and
Ivanovs (2021) for a review of recent developments in the literature of
multivariate extremes. Much of this line of research tries to connect
important ideas from modern statistics and machine learning to the context of
multivariate extremes. Our work falls in this category as we propose spectral
clustering as a tool for learning the dependence structure of multivariate
extremes.

Spectral clustering (Von Luxburg, 2007) and related techniques are very
popular and have found success in various applications such as parallel
computing (Hendrickson and Leland, 1995; Van Driessche and Roose, 1995), image
segmentation (Shi and Malik, 2000) and community detection (Rohe et al., 2011;
Lei and Rinaldo, 2015; Zhou and Amini, 2019). The central idea of spectral
clustering is to use the eigenvectors of the graph Laplacian matrix
constructed from an affinity graph between sample points in order to find
clusters in the data. Typically these are obtained by a $K$-means algorithm
that take these graph Laplacian eigenvectors as input. We follow this same
principle but use as input to our algorithm the angular parts of the
observations whose norms exceed a certain large threshold i.e., a standard
spectral clustering algorithm is applied to the graph built over the angular
parts of these extreme observations.

Because of the nature of the extreme events that we study, we leverage tools
from multivariate extreme value theory for analyzing the theoretical
properties of our spectral clustering algorithm. In particular, we use
multivariate regular variation as a modeling tool since it is closely
connected to asymptotic characterizations of multivariate extreme value
distributions (Resnick, 2007, 2018). While a precise definition of regular
variation is provided in Section 2, the basic idea is that a $d$-dimensional
random vector ${\bf X}$ is regularly varying if the distribution of the
angular part ${\bf X}/\|{\bf X}\|$ stabilizes (i.e., converges in
distribution) as the radial part $\|{\bf X}\|$ becomes large and that the
radial part has Pareto-like tails. The dependence structure is then governed
by the asymptotic distribution of the limiting angular part. In this paper, we
consider clustering of the angular parts, which live on a $d$-dimensional unit
sphere, of observations with large radii. Learning this measure is challenging
because of its multivariate nature and because only a small fraction of the
data is considered to be extremes, i.e., those observations whose radii are
sufficiently large, are retained for estimation. In contrast, standard
modeling approaches built on parametric models are hard to extend to larger
dimensions because of their lack of flexibility and computational complexity
Davison and Huser (2015).

We will explore the use of spectral clustering for learning the angular
measure. The performance of the algorithm critically depends on the properties
of the random graph that it takes as input. We will focus on $k$-nearest
neighbor graphs and hence a decision has to be made about the size of $k$ for
constructing the random graph. In this work we study this question by focusing
on a linear factor model. We characterize the asymptotic distribution of the
multivariate extremes generated from this model and show that their dependence
structure is captured by a discrete angular measure in the limit. We establish
a rate of convergence for the angular components of the extremes to their
discrete limits. This is a key step in deriving a theoretically valid range of
numbers of $k$-nearest neighbors for constructing a nearest neighbor graph
that one should consider in order to guarantee that spectral clustering can be
successfully used to learn the asymptotic angular measure.

From a methodological perspective, the work of Janßen and Wan (2020) is
perhaps the closest to our approach since they also provide a clustering
algorithm for extremes. Their method is however very different as it is based
on spherical $k$-means (Dhillon and Modha, 2001), a variant of $k$-means that
replaces the usual square loss minimization by an angular dissimilarity
measure minimization. The data-generating model we consider is a natural
factor model that can be viewed as a generalization of the max-linear model
considered in Janßen and Wan (2020). We characterize the limiting distribution
of the extremes in this model. We rigorously study the extremal nearest
neighbor graphs and show that their connected components can identify the
clusters of extremes of our factor model. By construction our algorithm is
computationally tractable and model agnostic, so it has a potential of working
well beyond the setting covered by our theory.

The rest of the paper is organized as follows. Section 2 provides some
background notions from multivariate regular variation necessary for our
analysis. Section 3 introduces the proposed spectral clustering algorithm for
extremes. In Section 4 we introduce our linear factor model (LFM) and derive
the asymptotic distribution of the angular components ${\bf X}/\|{\bf X}\|$ of
observations with high threshold exceedances i.e., observations ${\bf X}$ with
very large $\|{\bf X}\|$. In Section 5 we study the behavior of $k$-nearest
neighbor graphs constructed using a sample of extremes. Section 6 contains a
number of numerical examples that illustrate our proposed method. We show in
Section 6.1 that for a large range of values of $k$ the connected components
of the nearest neighbor graph consistently identify the clusters of extremes
arising from the linear factor model. This includes an examination of LFM with
added noise. The spectral clustering method is still able to estimate the
signal reasonably well. The good numerical performance of the method in the
LFM plus noise context suggests that it might work well in more general
settings. The spectral clustering method is also applied to an environmental
data set consisting of daily measurements of five air pollutants over both
winter and summer seasons. The analysis suggests that in modeling the
extremes, a LFM model with 5 clusters seems appropriate. Moreover, viewed as a
time series, the extremal dependence for O3 and NO2 does not extend beyond a
second-day time lag. Proofs of the technical results in the body of the paper
and their complements are contained in the appendix.

##  2\. Background on multivariate regular variation

Regular variation is often the starting point in modeling heavy-tailed data.
We will make regular use of this assumption throughout this work. A random
vector ${\bf X}=(X_{1},\ldots,X_{d})^{\top}$ is said to be regularly varying
with exponent $\alpha>0$ if for some norm $\|\cdot\|$ on $\mathbb{R}^{d}$ and
some probability measure $\Gamma$ on the unit sphere $\mathbb{S}^{d-1}$ in
$\mathbb{R}^{d}$, the following limits hold:

(2.1) $\displaystyle\lim_{r\to\infty}\mathbb{P}\bigl{(}{\bf X}/\|{\bf
X}\|\in\cdot\,\mid\,\|{\bf X}\|>r\bigr{)}\Rightarrow\Gamma(\cdot)$

and

(2.2) $\displaystyle\lim_{r\to\infty}\frac{\mathbb{P}\bigl{(}\|{\bf
X}\|>rx\bigr{)}}{\mathbb{P}\bigl{(}\|{\bf X}\|>r\bigr{)}}=x^{-\alpha}$

for all $x>0$, where $\Rightarrow$ denotes weak convergence on
$\mathbb{S}^{d-1}$. In other words, the law of the angular component ${\bf
X}/\|{\bf X}\|$ stabilizes as the radial component becomes large, and the
radial component is regularly varying (equation (2.2)) with index $\alpha$.
The limit probability measure $\Gamma$ is called the angular measure (or
spectral measure) and describes how likely the extremal observations are to
point in different directions. In other words, the angular measure describes
the limiting extremal angle for high threshold exceedances that correspond to
large $\|{\bf X}\|$. The support of this measure is particularly important
since it shows which directions of the extremes are feasible and which are not
feasible. Throughout the rest of paper we will take $\|\cdot\|$ to be the
Euclidean norm.

For example, if ${\bf X}$ has a spherically symmetric distribution and the
radius $\|{\bf X}\|$ has a Pareto distribution with index $\alpha$, then ${\bf
X}$ is regularly varying with angular measure that is uniform on
$\mathbb{S}^{d-1}$. In this case, the random vector is equally likely to have
extremes in any direction so we do not expect extremes to be clustered. On the
other hand, consider observations generated from a univariate MA(3) process
given by $Y_{t}=Z_{t}+.5Z_{t-1}-.6Z_{t-2}+1.5Z_{t-3}$, where $\\{Z_{t}\\}$ is
an iid sequence of symmetric stable random variables with index $\alpha=1.8$.
The bivariate vector ${\bf X}_{t}=(Y_{t},Y_{t-1})^{\top}$ is regularly varying
and the scatter plot of $Y_{t}$ vs $Y_{t-1}$ is displayed in the left panel of
Figure 1. Notice that for large values of $\|{\bf X}_{t}\|$, the points align
themselves on rays. In the right panel is a plot of ${\bf X}_{t}/\|{\bf
X}_{t}\|$ for those values of $\|{\bf X}_{t}\|$ that exceed the 99.8%
empirical quantile of the radii and are grouped in 10 clusters. In this
particular case, the spectral distribution consists of 10 point masses (5
pairs of symmetric point masses, indicated by arrows emanating from the
origin).

Figure 1. Scatter plot of $(Y_{t},Y_{t-1})$ for an MA(3) process (left);
spectral measure on $\mathbb{S}^{1}$ (right)

This simple example illustrates the challenge in finding meaningful low
dimensional regions supporting the extremes. In a series of papers (see Meyer
and Wintenberger (2019), Drees and Sabourin (2021), Cooley and Thibaud
(2019)), PCA-like analyses have been applied to finding low-dimensional
subspaces that contain the bulk of the support of the spectral measure. As
seen in this MA(3) example, such strategies might not be well suited for
extracting the key features in the extremes which do not necessarily live
neatly on a small number of subspaces. The approach taken here will be to use
spectral clustering to learn the angular measure. While machine learning ideas
provide guidance for thinking about and addressing multivariate extremes, the
very nature of rare events will require us to borrow ideas from the theory of
multivariate regular variation to analyze the extremal nearest neighbor graphs
used by our algorithm.

##  3\. Spectral clustering

In this section we describe how to construct random graphs based on a sample
of extremes and how to use such graphs to find clusters of extremes via a
simple spectral clustering algorithm.

###  3.1. Constructing random graphs

Starting with a sample of $d$-dimensional observations ${\bf
X}_{i},\,i=1,\ldots,n$, one first needs to identify the extremal part of the
sample, on which the extremal estimation will be performed. This is often done
by selecting a high threshold $u_{n}$ and assigning to the extremal part of
the sample the observations ${\bf X}_{i}$ satisfying $\|{\bf X}_{i}\|>u_{n}$.

Assume that $N_{n}$ observations ${\bf X}_{i}$ (with $i$ in some set
$\mathcal{V}_{n}$ of cardinality $N_{n}$) are in the extremal part of the
sample. Associated with each $i\in\mathcal{V}_{n}$, is the angular component
of the observation ${\bf X}_{i}/\|{\bf X}_{i}\|$ that lives on the unit
sphere. This allows us to think of the points in $\mathcal{V}_{n}$ as points
on the unit sphere, forming nodes in a simple graph. We connect nodes $i_{1}$
and $i_{2}$ by an edge according to a certain rule. One possible rule chooses
$\epsilon>0$ and connects $i_{1},i_{2}\in\mathcal{V}_{n}$ by an edge if

(3.1) $\displaystyle\rho\bigl{(}{\bf X}_{i_{1}}/\|{\bf X}_{i_{1}}\|,{\bf
X}_{i_{2}}/\|{\bf X}_{i_{2}}\|\bigr{)}\leq\epsilon\,.$

One often uses the usual Euclidean distance $\rho$ on the unit sphere
$\mathbb{S}^{d-1}$ in (3.1); but a distance function on the unit sphere could
also be used. The random set of edges $\mathcal{E}_{n}$ created in this
fashion define an $\epsilon$-neighborhood graph.

In what follows, we will focus on a different rule, leading to the $k$-Nearest
Neighbor graphs ($k$-NN graphs). This rule asserts that a node
$i_{1}\in\mathcal{V}_{n}$ is connected to a node $i_{2}\in\mathcal{V}_{n}$ if
the point on the unit sphere corresponding to $i_{2}$ is among the $k$-nearest
neighbors of the point corresponding to $i_{1}$, according to some distance
function. This definition leads to a directed graph because the neighborhood
relationship is not symmetric. There are two natural ways of making this graph
undirected. The first one is to connect $i_{1}$ and $i_{2}$ with an undirected
edge if either $i_{1}$ is among the $k$-nearest neighbors of $i_{2}$ or
$i_{2}$ is among the $k$-nearest neighbors of $i_{1}$. The second one connects
$i_{1}$ and $i_{2}$ only if both conditions are met, i.e., when $i_{1}$ and
$i_{2}$ are mutual nearest neighbors (hence the resulting graph is usually
called mutual $k$-nearest neighbor graph). Our main results apply to both
constructions. We work with weighted graphs, where we assign to the edges a
weight equal to the distance between the points on the unit sphere defining
the nodes. More specifically, we will take as input to our algorithm the
weighted adjacency matrix ${\bf
W}=[w_{i_{1}i_{2}}]_{i_{1},i_{2}\in\mathcal{V}_{n}}$ and

(3.2) $w_{i_{1}i_{2}}=\begin{cases}d\bigl{(}{\bf X}_{i_{1}}/\|{\bf
X}_{i_{1}},\|,{\bf X}_{i_{2}}/\|{\bf X}_{i_{2}}\|\bigr{)}&\mbox{if $i_{1}$ and
$i_{2}$ are connected,}\\\ 0,&\mbox{if $i_{1}$ and $i_{2}$ are not
connected.}\end{cases}$

When defining the weights in (3.2) $d$ is a certain kernel; a typical example
of such a kernel e.g., $d({\bf x},{\bf y})=\exp(-\|{\bf x}-{\bf y}\|)$, is
used in the examples of Section 6. In the following subsections, we describe
our algorithm and highlight the theoretical challenges.

###  3.2. The algorithm

The degree of a node $i\in\mathcal{V}_{n}$ is defined as

$d_{i}=\sum_{j\in\mathcal{V}_{n}}w_{ij}.$

The degree matrix $D$ is defined as the diagonal matrix with diagonal elements
$[d_{i}]_{i\in\mathcal{V}_{n}}$ and the normalized symmetric graph Laplacian
matrix is defined as

(3.3) $L=I-D^{-1/2}WD^{-1/2},$

where $I$ is the identity matrix. The spectral clustering algorithm of Ng et
al. (2002) proceeds as follows:

  1. (1)

Compute the first $m$ eigenvectors ${\bf u}_{1},\dots,{\bf u}_{m}$ of $L$
(i.e., the eigenvectors corresponding to the $m$ smallest eigenvalues of $L$)
and define an $N_{n}\times m$ matrix $U$ using these eigenvectors.

  2. (2)

Form an $N_{n}\times m$ matrix $V$ by normalizing the rows of $U$ to have unit
norm.

  3. (3)

Treating each of the $N_{n}$ rows of $V$ as a vector in $\mathbb{R}^{m}$,
cluster them into $m$ clusters $C_{1},\dots,C_{m}$ using the $K$-means
algorithm.

  4. (4)

Assign the original points ${\bf X}_{i}$ to cluster $C_{j}$ if and only if row
$i$ of the matrix $V$ was assigned to cluster $C_{j}$.

The motivation for this algorithm is described below.

###  3.3. Connected components, Laplacian and $k$-nearest neighbor graph.

We say that a subset $\mathcal{A}\subset\mathcal{V}_{n}$ of the vertices of a
graph is connected if any two vertices in $\mathcal{A}$ can be joined by a
path of edges such that all intermediate vertices also lie in $\mathcal{A}$.
If $\mathcal{A}$ is connected and there are no connections between
$\mathcal{A}$ and $\mathcal{V}_{n}\setminus\mathcal{A}$, then $\mathcal{A}$ is
called a connected component. It is well known that the number of connected
components of a graph $G$ is related to the spectrum of its symmetric graph
Laplacian. This is formalized in the following proposition (Von Luxburg, 2007,
Proposition 2).

######  Proposition 3.1.

Let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ be an undirected graph with non-
negative weights. Then the multiplicity $m$ of the eigenvalue $0$ of $L$
equals the number of connected components
$\mathcal{A}_{1},\dots,\mathcal{A}_{m}$ in the graph. The eigenspace of
eigenvalue $0$ is spanned by the indicator vectors
$\delta_{\mathcal{A}_{1}},\dots,\delta_{\mathcal{A}_{m}}$ of those components.

It follows from this result that if the spectral clustering algorithm is
applied to a graph with the number of connected components equal to the
parameter $m$ in the algorithm, then the algorithm will identify the connected
components. In Sections 4 and 5 we derive the asymptotic angular distribution
of multivariate extremes arising from a linear factor model and provide the
relevant asymptotic theory for the connected components of the $k_{n}$-nearest
neighbor graph constructed using the angular components of the extremes.
Specifically, it will be rigorously established that under certain conditions
the spectral clustering of the resulting graph consistently estimates the
support of the spectral measure of multivariate extremes arising from this
model.

##  4\.  Linear factor model and convergence of the angular components

We now introduce the generative model that we will be studying in this paper.
Let ${\bf X}$ be a $d$-dimensional random vector defined by the following
linear factor model (LFM)

(4.1) ${\bf X}=A{\bf Z}\,,$

where $A=[a_{ij}]_{i=1,\ldots d;j=1,\ldots p}$ is a $d\times p$ matrix of
nonnegative elements and ${\bf Z}$ is a $p$-dimensional random vector of
factors consisting of independent and identically distributed random
variables, that are either nonnegative or symmetric, and have asymptotically
Pareto tails, i.e.,

(4.2) $\mathbb{P}(Z_{1}>z)\sim cz^{-\alpha},\ \mbox{as $z\to\infty$}$

for some $\alpha>0$ and $c>0$. Note that we write $f(x)\sim g(x)$ as
$x\to\infty$ to mean that $\lim_{x\to\infty}f(x)/g(x)=1$. In the examples
section, we will add noise to the model in (4.1) in which case, the model
corresponds to a standard heavy-tailed linear factor model. One can think of
(4.1) as a linear version of the max-linear model studied in Janßen and Wan
(2020), which has the same spectral distribution. We will also relax the
assumption that the matrix $A$ is non-negative and will allow the noise to be
symmetric. In this case, the spectral distribution of the model is no longer
constrained to the positive quadrant of $\mathbb{S}^{d-1}$. Related max-linear
models have also been considered in the context of time series models for
extremes (Davis and Resnick, 1989; Hall et al., 2002) and more recently in the
context of structural equation models (Gissibl and Klüppelberg, 2018;
Klüppelberg and Lauritzen, 2019). The asymptotic Pareto assumption in (4.2)
can be weakened to regular variation, at least for the main results in this
section. However, this comes at the expense of assuming a more intrinsically
complex set of conditions on the choice of thresholding sequences. Additional
conditions, such as existence and properties of the density function of the
noise, are required for the proofs of the results in Section 5.

It follows immediately from (4.1) and (4.2) (see, for example, Basrak et al.
(2002), Proposition A.1) that ${\bf X}$ is a multivariate regularly varying
random vector satisfying (2.1); namely,

(4.3) $\lim_{x\to\infty}\mathbb{P}\left(\ \frac{{\bf X}}{\|{\bf
X}\|}\in\cdot~{}|~{}\|{\bf X}\|>x,\right)\Rightarrow\Gamma(\cdot)\,,$

where $\Rightarrow$ denotes weak convergence on the unit sphere
$\mathbb{S}^{d-1}$, $\Gamma$ is a discrete probability measure on
$\mathbb{S}^{d-1}$ that, in the nonnegative case, puts mass $\|{\bf
a}^{(k)}\|^{\alpha}/w$ at ${\bf a}^{(k)}/\|{\bf a}^{(k)}\|$ for
$k=1,\ldots,p$, where ${\bf a}^{(k)}=(a_{1k},a_{2k},\ldots,a_{dk})^{\top}\,,$
is the $k^{th}$ column of the matrix $A$ and

(4.4) $w=\sum_{k=1}^{p}\|{\bf a}^{(k)}\|^{\alpha}\,.$

In other words, $\Gamma$ has the representation

(4.5) $\Gamma(\cdot)=w^{-1}\sum_{k=1}^{p}\|{\bf
a}^{(k)}\|^{\alpha}\delta_{\frac{{\bf a}^{(k)}}{\|{\bf a}^{(k)}\|}}(\cdot)\,,$

where $\delta_{x}(\cdot)$ is the Dirac measure that puts unit mass at $x$. On
the other hand, in the symmetric case, $\Gamma$ puts mass $\|{\bf
a}^{(k)}\|^{\alpha}/2w$ at $\pm{\bf a}^{(k)}/\|{\bf a}^{(k)}\|$ for
$k=1,\ldots,p$. That is, the number of point masses in the symmetric case is
double of that number in the nonnegative case. 111Without much additional
effort, one could consider the case that the tails of $Z_{1}$ are balanced in
the sense that $\lim_{x\to\infty}\mathbb{P}(Z_{1}>x)/\mathbb{P}(|Z_{1}|>x)\to
p^{+}\in[0,1].$ The location of the point masses for $\Gamma$ would be exactly
the same as in the symmetric case, but with mass $p^{\pm}\|{\bf
a}^{(k)}\|^{\alpha}/w$ at ${\bf c}_{k,\pm}$ defined in (4.7), where
$p^{-}=1-p^{+}$.

Based on a random sample of iid copies of ${\bf X}_{1},\ldots,{\bf X}_{n}$ of
${\bf X}$ as above, we construct an estimate of the location of the point
masses that comprise $\Gamma$, i.e.,

(4.6) ${\bf c}_{k}=\frac{{\bf a}^{(k)}}{\|{\bf a}^{(k)}\|}\,,~{}k=1,\ldots,p$

in the nonnegative case, and

(4.7) ${\bf c}_{k,\pm}=\frac{\pm{\bf a}^{(k)}}{\|{\bf
a}^{(k)}\|}\,,~{}k=1,\ldots,p$

in the symmetric case. Note that these ${\bf c}_{k}$ are not necessarily
distinct. Intuitively, for large $n$, the angular parts ${\bf X}_{i}/\|{\bf
X}_{i}\|$ of the sample for which $\|{\bf X}_{i}\|$ is large, will cluster
around these ${\bf c}_{k}$. In fact, we formalize this intuition and provide a
rate of convergence for the limiting extremal angles with high threshold
exceedances in the next theorem. This will be a key ingredient in our
convergence analysis of extremal $k$-NN graphs. For the ease of notation we
will prove the following result in the nonnegative case; the symmetric case
follows by simply doubling the number of points on the sphere.

######  Theorem 4.1.

If $(u_{n})$ is a sequence converging to infinity as $n\to\infty$, then, in
the nonnegative case, for any $j=1,\ldots,p$, the conditional law of

$u_{n}\bigl{(}{\bf X}/\|{\bf X}\|-{\bf c}_{j}\bigr{)}$

given $\|{\bf X}\|>u_{n},\,Z_{j}>u_{n}/w^{1/\alpha}$, ($w$ defined in (4.4))
converges weakly to the law of

$\frac{1}{\|{\bf
a}^{(j)}\|^{2}W_{\alpha}}\left(S^{*}_{1,-j},\ldots,S^{*}_{d,-j}\right)^{\top},\,$

where $W_{\alpha}$ has a Pareto distribution (i.e.
$\mathbb{P}(W_{\alpha}>x)=x^{-\alpha},\,x\geq 1$) that is independent of
$Z_{1},\ldots,Z_{p}$,

(4.8)
$S^{*}_{l,-j}=\sum_{i=1}^{d}\left(a_{ij}^{2}X_{l,-j}-a_{lj}a_{ij}X_{i,-j}\right)\,,$

and

(4.9) $X_{i,-j}=X_{i}-a_{ij}Z_{j}=\sum_{\begin{subarray}{c}m=1,\\\ m\neq
j\end{subarray}}^{p}a_{im}Z_{m}\,.$

###### Proof.

We start by observing that the conditional law of

(4.10)
$\bigl{(}Z_{1},\ldots,Z_{j-1},Z_{j}/u_{n},Z_{j+1},\ldots,Z_{p}\bigr{)}^{\top}$

given $\|{\bf X}\|^{2}>u_{n}^{2}$, $Z_{j}>u_{n}/w^{1/\alpha}$, converges in
distribution, as $n\to\infty$, to the law of

(4.11)
$\bigl{(}Z_{1},\ldots,Z_{j-1},W_{\alpha}/w_{j},Z_{j+1},\ldots,Z_{p}\bigr{)}^{\top},$

where $w_{j}=\|{\bf a}^{(j)}\|$ The main ingredients in establishing this
result is to note that $Z_{j}^{2}$ is regularly varying with index $\alpha/2$
while for $i\neq j$, $Z_{i}Z_{j}$ is regularly varying with index $\alpha$
(see Embrechts and Goldie (1980); Theorem 3). Moreover, from the convolution
closure property for sums of independent regularly varying random variables,
it follows easily that

$\displaystyle\mathbb{P}(Z_{j}>u_{n}x\,,\sum_{k=1}^{d}(\sum_{i=1}^{p}a_{ki}Z_{i})^{2}>u_{n}^{2},\,)$
$\displaystyle\sim$
$\displaystyle\mathbb{P}(Z_{j}>u_{n}x,\sum_{i=1}^{p}w_{i}^{2}Z_{i}^{2}>u_{n}^{2})$
$\displaystyle\sim$
$\displaystyle\sum_{i=1}^{p}\mathbb{P}(Z_{j}>u_{n}x,w_{i}^{2}Z_{i}^{2}>u_{n}^{2})$
$\displaystyle\sim$ $\displaystyle\mathbb{P}(Z_{j}>u_{n}x)\,.$

Now to finish the proof of (4.11), we use these relations and note that for
$x\geq 1/w_{j}$ (and hence $x\geq 1/w^{1/\alpha}$),

$\displaystyle\mathbb{P}(Z_{j}>u_{n}x\,|\,\|{\bf
X}\|^{2}>u_{n}^{2},\,Z_{j}^{2}>u_{n}^{2}/w^{2/\alpha})$ $\displaystyle\sim$
$\displaystyle\frac{\mathbb{P}(Z_{j}>u_{n}x,\sum_{i=1}^{p}w_{i}^{2}Z_{i}^{2}>u_{n}^{2},Z_{j}^{2}>u_{n}^{2}/w^{2/\alpha})}{\mathbb{P}(\sum_{i=1}^{p}w_{i}^{2}Z_{i}^{2}>u_{n}^{2},Z_{j}^{2}>u_{n}^{2}/w^{2/\alpha})}$
$\displaystyle\sim$
$\displaystyle\frac{\mathbb{P}(Z_{j}>u_{n}x)}{\mathbb{P}(Z_{j}>u_{n}/w_{j})}$
$\displaystyle\to$ $\displaystyle w_{j}^{-\alpha}x^{-\alpha}$ $\displaystyle=$
$\displaystyle\mathbb{P}(W_{\alpha}>w_{j}x)\,.$

We have

$\displaystyle u_{n}\bigl{(}{\bf X}/\|{\bf X}\|-{\bf c}_{j}\bigr{)}$
$\displaystyle=u_{n}\left(\frac{\left(\sum_{m=1}^{p}a_{1m}Z_{m},\ldots,\sum_{m=1}^{p}a_{dm}Z_{m}\right)^{\top}}{\|{\bf
X}\|}-\frac{{\bf a}^{(j)}}{\|{\bf a}^{(j)}\|}\right)$
$\displaystyle=(V_{1},\ldots,V_{d})^{\top},$

say. For $l=1,\ldots,d$,

$\displaystyle V_{l}$
$\displaystyle=u_{n}\frac{w_{j}\sum_{m=1}^{p}a_{lm}Z_{m}-a_{lj}\|{\bf
X}\|}{w_{j}\|{\bf X}\|}$ (4.12)
$\displaystyle=u_{n}\frac{w_{j}^{2}\left(\sum_{m=1}^{p}a_{lm}Z_{m}\right)^{2}-a_{lj}^{2}\|{\bf
X}\|^{2}}{w_{j}\|{\bf X}\|\left(w_{j}\sum_{m=1}^{p}a_{lm}Z_{m}+a_{lj}\|{\bf
X}\|\right)}\,.$

We note that the numerator of (4.12) reduces to the following expression where
the $Z_{j}^{2}$ terms cancel out

(4.13) $\displaystyle Num_{l}:$ $\displaystyle=$ $\displaystyle
2a_{lj}Z_{j}\left(\|{\bf
a}^{(j)}\|^{2}X_{l,-j}-a_{lj}\sum_{i=1}^{d}a_{i,j}X_{i,-j}\right)+\|{\bf
a}^{(i)}\|^{2}X_{l,-j}^{2}-a_{lj}^{2}\sum_{i=1}^{d}X_{i,-j}^{2}$
$\displaystyle=$ $\displaystyle
2a_{lj}Z_{j}\left(\sum_{i=1}^{d}\left(a_{ij}^{2}X_{l,-j}-a_{lj}a_{ij}X_{i,-j}\right)\right)+\|{\bf
a}^{(i)}\|^{2}X_{l,-j}^{2}-a_{lj}^{2}\sum_{i=1}^{d}X_{i,-j}^{2}\,,$

where $X_{i,-j}$ is as defined in (4.9). The denominator of (4.12) is handled
in a similar way, but this time the $Z_{j}^{2}$ terms do not cancel. Indeed,
since

$\displaystyle\|{\bf X}\|$
$\displaystyle=\sqrt{\sum_{k=1}^{d}(a_{kj}^{2}Z_{j}^{2}+2a_{kj}Z_{j}X_{k,-j}+X_{k,-j}^{2})}$
$\displaystyle=\sqrt{\sum_{k=1}^{d}a_{kj}^{2}Z_{j}^{2}}\sqrt{1+\frac{\sum_{k=1}^{d}(X_{k,-j}^{2}+2a_{kj}Z_{j}X_{k,-j})}{\sum_{k=1}^{d}a_{kj}^{2}Z_{j}^{2}}}$
$\displaystyle=w_{j}|Z_{j}|\sqrt{1+\frac{\sum_{k=1}^{d}(X_{k,-j}^{2}+2a_{kj}Z_{j}X_{k,-j})}{w_{j}^{2}Z_{j}^{2}}},$

we can write

$Den_{l}=w_{j}^{2}|Z_{j}|(1+o_{p}(1))\left(w_{j}a_{lj}Z_{j}+R+a_{lj}w_{j}|Z_{j}|(1+o_{p}(1))\right),$

where $R$ is a linear function in the variables $Z_{1},\ldots,Z_{p}$ in which
$Z_{j}$ does not appear, and $o_{p}(1)$ goes to zero in probability given
$Z_{j}>u_{n}/w^{1/\alpha}$. We view

$V_{l}=u_{n}\frac{Num_{l}}{Den_{l}}=\frac{Num_{l}/u_{n}}{Den_{l}/u_{n}^{2}}$

as a ratio of two continuous real-valued functions of the random vector in
(4.10) (plus a vanishing term in $Dem_{l}$), so that the random vector
$(V_{1},\ldots,V_{d})^{\top}$ becomes a $d$-dimensional vector of such ratios.
By the continuous mapping theorem the random vector
$(V_{1},\ldots,V_{d})^{\top}$ converges weakly to the $d$-dimensional vector
of the ratios of the corresponding functions applied to the random vector in
(4.11). These result in

$2a_{lj}(W_{\alpha}/w_{j})S^{*}_{l,-j}$

in the case of $Num_{l}$ and in

$W_{\alpha}^{2}\,2a_{lj}w_{j}$

in the case of $Dem_{l}$. Putting everything together produces the claim. ∎

As explained above, the following corollary is an immediate consequence of
Theorem 4.1.

######  Corollary 4.2.

Let a sequence of levels $(u_{n})$ converging to infinity. Then, in the
symmetric case, in the notation of Theorem 4.1, for any $j=1,\ldots,p$, the
conditional law of

$u_{n}\bigl{(}{\bf X}/\|{\bf X}\|-\pm{\bf c}_{j}\bigr{)}$

given $\|{\bf X}\|>u_{n},\,\pm Z_{j}>u_{n}/w^{1/\alpha}$, converges weakly to
the law of

$\frac{\pm 1}{\|{\bf
a}^{(j)}\|^{2}W_{\alpha}}\left(S^{*}_{1,-j},\ldots,S^{*}_{d,-j}\right)^{\top}.$

######  Remark 4.3.

In the nonnegative case, Theorem 4.1 addresses the conditional convergence of
${\bf X}/\|{\bf X}\|$, given $\|{\bf X}\|>u_{n}$, $Z_{j}>u_{n}/w^{1/\alpha}$,
to the location ${\bf c}_{j}$ of the corresponding atom of the spectral
measure. It is also possible to address a conditional convergence to the mass
$w^{-1}\|{\bf a}^{(j)}\|^{\alpha}$ of this atom. Indeed,

(4.14) $\mathbb{P}\bigl{(}Z_{j}>u_{n}/w^{1/\alpha}\big{|}\|{\bf
X}\|>u_{n}\bigr{)}\to w^{-1}\|{\bf a}^{(j)}\|^{\alpha}$

as $n\to\infty$. To see this, write

$\displaystyle\mathbb{P}\bigl{(}Z_{j}>u_{n}/w^{1/\alpha}\big{|}\|{\bf
X}\|>u_{n}\bigr{)}=\frac{u_{n}^{\alpha}\mathbb{P}\bigl{(}Z_{j}>u_{n}/w^{1/\alpha},\|{\bf
X}\|>u_{n}\bigr{)}}{u_{n}^{\alpha}\mathbb{P}\bigl{(}\|{\bf
X}\|>u_{n}\bigr{)}},$

and the numerator converges to $c\|{\bf a}^{(j)}\|^{\alpha}$, while the
denominator converges to $cw$. If one strengthens the asymptotic Pareto tails
assumption (4.2) to include the rate of convergence to the limit, then one
would able to establish the rate of convergence in (4.14) as well. The
situation is similar in the symmetric case. We do not pursue this in the
present paper.

We now explore the connection between large values of the underlying factors
$Z_{i1},\dots,Z_{ip}$ and large values of $\|{\bf X}_{i}\|$. We will see that
under certain conditions, high threshold exceedances of $\|{\bf X}_{i}\|$ are
generated by only one underlying factor $Z_{ij}$, $j=1,\dots,p.$ This will be
important for our analysis of extremal $k$-NN graphs which will require
additional assumptions on the rate of growth of $u_{n}$. Since
$\frac{\alpha+2}{\alpha(\alpha+3)}<\alpha^{-1}$, we can further impose that
the sequence $(u_{n})$ satisfies the growth conditions

(4.15) $n^{-1/\alpha}u_{n}\to
0~{}~{}\mbox{and}~{}~{}n^{-(\alpha+2)/(\alpha(\alpha+3))}u_{n}\to\infty\,,$

as $n\to\infty$. Also note that we may choose a further sequence $(h_{n})$
such that

(4.16)
$h_{n}\to\infty,\,h_{n}=o(u_{n}),\,h_{n}=o\bigl{(}u_{n}^{(\alpha+1)/2}n^{-1/2}\bigr{)},\,n^{-1/\alpha}u_{n}h_{n}\to\infty$

as $n\to\infty$. Indeed, the choice

$h_{n}=u_{n}^{(\alpha-1)/4}n^{(2-\alpha)/(4\alpha)}$

works for this purpose.

For $n=1,2,\ldots,$ we define the set of indexes corresponding to extreme
observations

(4.17) $\mathcal{I}_{n}=\bigl{\\{}i=1,\ldots,n:\,\|{\bf
X}_{i}\|>u_{n}\bigr{\\}},$

and denote its cardinality by $N_{n}=$card$(\mathcal{I}_{n})$. From (4.2) and
(4.15), we see that the mean and variance of $N_{n}/(nu_{n}^{-\alpha})$
converge to $cw$ and $0$, respectively and hence that

(4.18) $N_{n}/(nu_{n}^{-\alpha})\stackrel{{\scriptstyle
P}}{{\to}}cw,\mbox{~{}as $n\to\infty$.}$

The following lemma connects exceedances of $u_{n}$ by $\|{\bf X}_{i}\|$ with
exceedances of $h_{n}$ by $Z_{ij},\,j=1,\ldots,p$.

######  Lemma 4.4.

Let $(h_{n})$ be a sequence satisfying (4.16) and consider the event

${B}_{n}=\bigl{\\{}\text{for any $i\in\mathcal{I}_{n}$ at most one of
$Z_{im},\,m=1,\ldots,p$ exceeds $h_{n}$}\bigr{\\}}.$

Then $\mathbb{P}({B}_{n})\to 1$ as $n\to\infty$.

###### Proof.

Note that (4.1) implies that the $m^{th}$ component of the $i^{th}$
observation is of the form

(4.19) $\displaystyle
X_{im}=\sum_{j=1}^{p}a_{mj}Z_{ij},\,m=1,\ldots,d;~{}i=1,\ldots,n\,,$

where $Z_{i1},\dots,Z_{ip}$ are iid random variables with asymptotic Pareto
tails (4.2). Denote

(4.20) $a^{*}=d^{1/2}\max\\{a_{mj},\,{m=1,\ldots,d;\,j=1,\ldots,p}\\}$

Since $u_{n}>h_{n}$ for $n$ large, we have

$\displaystyle\mathbb{P}({B}_{n}^{c})\leq\sum_{i=1}^{n}\mathbb{P}\left(\sum_{k=1}^{d}\bigl{(}X_{ik}\bigr{)}^{2}>u_{n}^{2},\,Z_{im}>h_{n}\
\text{for two or more of }\ m=1,\ldots,p\right)$ $\displaystyle\leq
n\mathbb{P}\left(a^{*}\max_{k=1,\ldots,p}Z_{1k}>u_{n},\ Z_{1m}>h_{n}\
\text{for two or more of }\ m=1,\ldots,p\right)$
$\displaystyle\leq\sum_{k=2}^{p}{p\choose
k}n\mathbb{P}(Z_{1}>h_{n})\mathbb{P}\left(Z_{1}>u_{n}/a^{*}\right)\to 0$

by the last property in (4.16) and (4.2). This proves the lemma.

∎

Equipped with Lemma 4.4 we can now proceed to bound the distance between the
observed angular parts of the multivariate extremes and their corresponding
theoretical asymptotic atoms. Assume, for a moment, that we are in the
nonnegative case. We already know that for large $n$, we have that for every
$i\in\mathcal{I}_{n}$ one of the values of $Z_{im},\,m=1,\ldots,p$ must exceed
$u_{n}/a^{*}$ and all other values of these variables cannot exceed $h_{n}$.
We now define the sets of indexes corresponding to extremes generated by each
of the individual factors i.e. we define for $j=1,\ldots,p$

(4.21) $\mathcal{I}_{n}^{(j)}=\left\\{i=1,\ldots,n:\,\|{\bf
X}_{i}\|>u_{n},\,Z_{ij}>u_{n}/a^{*}\right\\}.$

Consequently,

(4.22) $\mathcal{I}_{n}=\bigcup_{j=1}^{p}\mathcal{I}_{n}^{(j)}$

and by Lemma 4.4 for large $n$ this is a disjoint union with probability
tending to one. Let $N_{n}^{(j)}$ be the cardinality of
$\mathcal{I}_{n}^{(j)}$, $j=1,\ldots,p$. Using the fact that $a^{*}\geq\|{\bf
a}^{(j)}\|$ for $j=1,\ldots,p$, the same argument as in (4.18) shows that
$j=1,\ldots p$,

(4.23) $N_{n}^{(j)}/(nu_{n}^{-\alpha})\stackrel{{\scriptstyle
P}}{{\to}}c\|{\bf a}^{(j)}\|^{\alpha},\mbox{~{}as $n\to\infty$.}$

We enumerate ${\bf X}_{i}/\|{\bf X}_{i}\|,\,i\in\mathcal{I}_{n}$ as ${\bf
Y}_{i},\,i=1,\ldots,N_{n}$, a sample on $\mathbb{S}^{d-1}$ of random size
$N_{n}$. For each $j=1,\ldots,p$, we enumerate ${\bf X}_{i}/\|{\bf
X}_{i}\|,\,i\in\mathcal{I}_{n}^{(j)}$ as ${\bf
Y}_{i}^{(j)},\,i=1,\ldots,N_{n}^{(j)}$, a sample on $\mathbb{S}^{d-1}$ of
random size $N_{n}^{(j)}$. It is straightforward (if a bit tedious) to check
the following result.

######  Lemma 4.5.

For large $n$, on the event ${B}_{n}$, for $i=1,\ldots,N_{n}^{(j)}$,

(4.24) $\displaystyle\Big{\|}{\bf Y}_{i}^{(j)}-{\bf
c}_{j}\Big{\|}\leq\frac{8(a^{*})^{2}}{\|{\bf
a}^{(j)}\|^{\alpha}}\frac{h_{n}}{u_{n}},\ \ j=1,\ldots,p,$

where the ${\bf c}_{j}$ are as defined in (4.7).

The situation in the symmetric case is, of course, completely analogous. It
follows from the definition of $h_{n}$ in (4.16) and (4.24) that the angular
components of the extremes are clustered around the centers ${\bf c}_{j}$. The
results in the next section build on Lemma 4.5 and provide sufficient
conditions for our extremal spectral clustering algorithm to be consistent.
For this, we provide a careful asymptotic analysis of the extremal $k$-NN
graph used by the algorithm.

##  5\.  Asymptotic analysis of the connected components of the extremal
$k$-NN graph

Our analysis consists of two main components. The first one is to show that
the extremes generated by different factors will belong to different
components of the extremal $k$-NN graph as long as the cluster centers
corresponding to the underlying factors are different. The second part will be
to argue that all the extremes generated by an underlying factor will also
belong to the same component of the extremal $k$-NN graph. This second step
turns out to be the more technical one in our analysis and we will only
establish this result for $d=2$. Along the way we derive a few intermediate
results that we also highlight in order to better explain the key ingredients
of our argument. Going forward, in our proofs, $c>0$ represents a finite and
non-zero constant whose value may change from line-to-line. In the sequel we
will assume, without further comments, that the sequence $(u_{n})$ satisfies
(4.15). The first step of our program is covered by the following proposition.

######  Proposition 5.1.

Suppose that $k_{n}=o\bigl{(}nu_{n}^{-\alpha}\bigr{)}$ as $n\to\infty$. Then
there is a sequence $({B}_{n,1})$ of events with $\mathbb{P}({B}_{n,1})\to 1$
as $n\to\infty$ such that, for all $n$ large enough, on the event ${B}_{n,1}$,
any two points ${\bf Y}_{i_{1}}^{(j_{1})}$ and ${\bf Y}_{i_{2}}^{(j_{2})}$,
$i_{1}=1,\ldots,N_{n}^{(j_{1})}$, $i_{2}=1,\ldots,N_{n}^{(j_{2})}$ will belong
to two different connected components of the $k_{n}$-NN graph if ${\bf
c}_{j_{1}}\not={\bf c}_{j_{2}}$.

###### Proof.

Define

(5.1) ${B}_{n,1}={B}_{n}\cap\bigl{\\{}N_{n}^{(j)}>k_{n}\ \ \text{for
$j=1,\ldots,p$}\bigr{\\}}.$

By Lemma 4.4, (4.23) and the assumption on $k_{n}$ , $\mathbb{P}({B}_{n,1})\to
1$ as $n\to\infty$. By (4.24) and the triangle inequality, on the events
${B}_{n,1}$, any point ${\bf Y}_{i_{1}}^{(j_{1})}$ has at least $k_{n}$
neighbours of the type ${\bf Y}_{i}^{(j_{1})}$,
$i=1,\ldots,N_{n}^{(j_{1})},\,i\not=i_{1}$, that are within distance of
$c\cdot h_{n}/u_{n}$ from it. On the other hand, by (4.24) and the triangle
inequality, its distance from any point ${\bf Y}_{i_{2}}^{(j_{2})}$ with ${\bf
c}_{j_{1}}\not={\bf c}_{j_{2}}$ cannot be smaller than

$\|{\bf c}_{j_{1}}-{\bf c}_{j_{2}}\|-c\cdot h_{n}/u_{n}.$

Therefore, for large $n$, the latter point cannot be among the $k_{n}$-nearest
neighbours of ${\bf Y}_{i_{1}}^{(j_{1})}$. ∎

We now embark on the second step of our program and establish that, at least
in the case $d=2$, under appropriate conditions, the points ${\bf
Y}_{i}^{(j)}$, $i=1,\ldots,N_{n}^{(j)}$, belong, with high probability, to the
same connected component in the $k_{n}$-NN graph. We start by investigating
the deviations of these points from the center of the cluster, ${\bf c}_{j}$,
defined in (4.7). Since the points ${\bf Y}_{i}^{(j)}$,
$i=1,\ldots,N_{n}^{(j)}$ are treated as independent, the following result is
essentially immediate from Theorem 4.1.

######  Lemma 5.2.

In the nonnegative case, for any $j=1,\ldots,p$, the conditional law of

$u_{n}\bigl{(}{\bf Y}_{1}^{(j)}-{\bf c}_{j}\bigr{)}$

given $N_{n}^{(j)}\geq 1$, converges weakly to the law of

$\frac{1}{\|{\bf
a}^{(j)}\|^{2}W_{\alpha}}\left(S^{*}_{1,-j},\ldots,S^{*}_{d,-j}\right)^{T}$

that is specified in the statement of Theorem 4.1. An analogous statement
holds in the symmetric case.

######  Remark 5.3.

It is a straightforward calculation to check that, if $j=1,\ldots,d$,

(5.2)
$\sum_{l=1}^{d}a_{lj}S_{l,-j}^{*}=\sum_{l=1}^{d}\sum_{i=1}^{d}\left(a_{lj}a_{ij}^{2}X_{l,-j}-a_{lj}^{2}a_{ij}X_{i,-j}\right)=0\,.$

Therefore, the normalized deviations of the points ${\bf Y}_{i}^{(j)}$,
$i=1,\ldots,N_{n}^{(j)}$ from the center of the $j$th cluster are, in the
limit, supported by a $(d-1)$-dimensional subspace.

Using the information in Lemma 5.2 we now proceed to prove that under
appropriate conditions, the points ${\bf Y}_{i}^{(j)}$,
$i=1,\ldots,N_{n}^{(j)}$, belong, with high probability, to the same connected
component of the $k_{n}$-NN graph. We will need some additional notation in
order to state the result. For a fixed $j=1,\ldots,d$ we write for
$m=1,\ldots,d$,

(5.3) $\displaystyle
Y_{im}^{(j)}=\sum_{l=1}^{p}a_{ml}Z^{(*,j)}_{il},\,i=1,\ldots,N_{n}^{(j)};$

the notation should be compared with (4.19). That is,
$\\{(Z^{(*,j)}_{i1},\ldots,Z^{(*,j)}_{ip})^{\top},i=1,\ldots,\\}$ are iid
random vectors distributed according to the conditional distribution of
$(Z_{1},\ldots,Z_{p})^{\top}$ given $\|{\bf
X}\|>u_{n},Z_{j}>u_{n}/w^{1/\alpha}$. Since the connectivity of any nearest
neighbor graph is not affected by shifting and scaling, it is sufficient to
consider the $k_{n}$-NN graph constructed on the deviations of the points
${\bf Y}_{i}^{(j)}$, $i=1,\ldots,N_{n}^{(j)}$ from the cluster center.

Continuing with the notation used in the proof of Theorem 4.1 we isolate the
main term in the deviations from the cluster center by writing

(5.4) $\displaystyle u_{n}\bigl{(}{\bf Y}_{i}^{(j)}-{\bf c}_{j}\bigr{)}$
$\displaystyle=\left(S_{1,-j}^{(*,j)},\ldots,S_{d,-j}^{(*,j)}\right)^{\top}/\bigl{(}w_{j}^{2}(Z^{(*,j)}_{ji}/u_{n})\bigr{)}$
$\displaystyle\quad\quad+\left[u_{n}\bigl{(}{\bf Y}_{i}^{(j)}-{\bf
c}_{j}\bigr{)}-\left(S^{(*,j)}_{1,-j},\ldots,S^{(*,j)}_{d,-j}\right)^{\top}/\bigl{(}w_{j}^{2}(Z^{(*,j)}_{ji}/u_{n})\bigr{)}\right]$
$\displaystyle={\bf M}^{(i)}+{\bf D}^{(i)},\quad i=1,\ldots,N_{n}^{(j)},$

where
$S_{l,-j}^{(*,j)}=Y_{il}^{(j)}-a_{lj}Z_{ij}^{(*,j)}=\sum_{\begin{subarray}{c}m=1,\\\
m\neq j\end{subarray}}^{p}a_{lm}Z_{im}^{(*,j)}\,$ is analogous to (4.9). In
the case $d=2$, it follows from (5.2) that for some nonzero deterministic
vector ${\bf b}$ in $\mathbb{R}^{2}$,

${\bf
M}^{(i)}=\frac{1}{w_{j}^{2}}\frac{S_{2,-j}^{(*,j)}}{Z^{(*,j)}_{ij}/u_{n}}{\bf
b},\ i=1,\ldots,N_{n}^{(j)}.$

For notational simplicity we continue the discussion with $j=1$, and in this
case these are essentially univariate iid random variables with the
distribution of

(5.5) $T_{n}=\frac{a_{22}Z_{2}+\cdots+a_{2p}Z_{p}}{w_{1}^{2}Z_{1}/u_{n}}$

given

(5.6)
$(a_{11}Z_{1}+\cdots+a_{1,p}Z_{p})^{2}+(a_{21}Z_{1}+\cdots+a_{2p}Z_{p})^{2}>u_{n}^{2},\
\ Z_{1}>u_{n}/w^{1/\alpha}\,.$

Finally, we let $F_{T_{n}}$ denote the conditional law of $T_{n}$ in (5.5)
given the conditions in (5.6). For technical reasons we require further
conditions on the latent factors in our subsequent results. We assume that the
generic noise variable $Z$ in (4.1) and (4.2) is positive or symmetric, and
has a probability density function $f_{Z}$ such that

(5.7) $f_{Z}(z)\ \ \text{is bounded away from 0 on compact intervals and
bounded from above,}$

and

(5.8) $B^{-1}z^{-(\alpha+1)}\leq f_{Z}(z)\leq Bz^{-(\alpha+1)},$

$\alpha>1$, for all $z\geq z_{0}$, some $B\geq 1$.

The following lemma shows that the conditional density function
$f_{T_{n}}(t)=\partial F_{T_{n}}(t)/\partial t$ enjoys some useful regularity
properties.

######  Lemma 5.4.

Assume (5.7) and (5.8). For $\alpha>1$ the density function $f_{T_{n}}$ is
such that:

  1. (i)

There exists an $G\in(0,\infty)$ such that for all $n$ large enough,
$f_{T_{n}}(t)\leq G$ for all $t$.

  2. (ii)

$f_{T_{n}}$ is uniformly bounded from below on compact intervals, uniformly in
$n$.

  3. (iii)

There is a constant $D\geq 1$ and a number $t_{0}\geq 0$ such that
$D^{-1}t^{-(\alpha+1)}\leq f_{T_{n}}(t)\leq Dt^{-(\alpha+1)}$ uniformly for
all $n$ large enough and all $t\geq t_{0}$.

It is clear that an analogous result holds for the appropriate conditional
densities in the symmetric case.

The following intermediate result is the key ingredient for completing our
analysis of the connected component of the extremal $k_{n}$-NN graph, at least
in the case $d=2$, assuming certain regularity conditions on the noise
variables.

######  Lemma 5.5.

Assume (5.7), (5.8) and let $d=2$, $\tau>1$ and consider the random variable
$m_{n}$ defined by

(5.9) $m_{n}=N_{n}^{(1)}/\lceil\tau\log N_{n}^{(1)}\rceil,\ \text{so that by
\eqref{e:In.size}}\ \
m_{n}\sim\frac{cw_{1}nu_{n}^{-\alpha}}{\tau\log(nu_{n}^{-\alpha})},\ \
n\to\infty.$

Define the intervals

(5.10)
$I_{i,n}=\left(F_{T_{n}}^{-1}\bigl{(}(i-1)/m_{n}\bigr{)},F_{T_{n}}^{-1}\bigl{(}i/m_{n}\bigr{)}\right),\
i=1,\ldots,m_{n}$

as well as intervals along vector ${\bf b}$ by

$J_{i,n}=I_{i,n}{\bf b},\ i=1,\ldots,N_{n}^{(1)}\,.$

Then, on an event with probability tending to one, there is a finite number
$K_{0}=1,2,\ldots$ such that for all $n$ large enough and all
$i=2,\ldots,m_{n}-I_{0}$, every point in $J_{i,n}$ is closer to every point in
the intervals $J_{i-1,n}$ and $J_{i+1,n}$ than to any point in an interval
$J_{i_{1},n}$ with $|i-i_{1}|>K_{0}$.

The proofs of these two lemmas are contained in the Appendix.

We are now ready to state the main result showing that the extremes generated
from the same underlying factor will also belong to the same connected
component of the extremal $k_{n}$-NN graph with probability tending to one,
under appropriate regularity conditions. This time, for simplicity, we only
consider the symmetric case.

######  Theorem 5.6.

Assume (5.7), (5.8) and let $d=2$. Then, if $k_{n}>{G}\log n$ with large
enough ${G}>0$, there is a sequence $({B}_{n,2})$ of events with
$\mathbb{P}({B}_{n,2})\to 1$ as $n\to\infty$ such that, for all $n$ large
enough, on the event ${B}_{n,2}$, any two points ${\bf Y}_{i_{1}}^{(j)}$ and
${\bf Y}_{i_{2}}^{(j)}$, $i_{1}=1,\ldots,N_{n}^{(j)}$,
$i_{2}=1,\ldots,N_{n}^{(j)}$ will belong to the same connected component of
the $k_{n}$-NN graph for any $j=1,\ldots,p$.

The proof of the theorem has been relegated to the appendix.

It follows from Proposition 5.1 and Theorem 5.6 that, with probability tending
to one as $n\to\infty$, the extremal $k_{n}$-NN graph obtained from a sample
drawn from (4.1) will have exactly $m\leq p$ connected components
corresponding to the $m$ distinct asymptotic point masses (4.7) of the model.
In other words, the extremal $k_{n}$-NN graph consistently identifies the
underlying clusters through its connected components. This in turn implies
that spectral clustering will be consistent by Proposition 3.1. We have
therefore shown the following main practical result.

######  Corollary 5.7.

Assume (5.7), (5.8), $d=2$, $k_{n}=o(nu_{n}^{-\alpha})$ and $k_{n}>G\log n$.
Then, spectral clustering will consistently identify the clusters of extremes
arising from the linear factor model.

######  Remark 5.8.

In practice consistent clustering can be achieved by taking $k_{n}>G_{0}\log
N_{n}$ for some $G_{0}>0$ and $k_{n}=o(N_{n})$ since
$N_{n}/(nu_{n}^{-\alpha})\stackrel{{\scriptstyle P}}{{\to}}cw$, as
$n\to\infty$. In our experiments we chose $k_{n}=\lceil\frac{N_{n}}{\tau\log
N_{n}}\rceil+1$ for some $\tau>1$.

Corollary 5.7 suggests a simple strategy for estimating the asymptotic angular
measure of the extremes generated from the linear factor model (4.1). Assume
we run spectral clustering on the extremal $k_{n}$-NN graph. Then we can
denote by $\hat{\mathcal{I}}_{n}^{(j)}$ the set of indices corresponding to
the $j$th cluster found by the algorithm for $j=1,\dots,m$. With these sets we
can define $\hat{N}_{n}^{(j)}=\mbox{card}(\mathcal{I}_{n}^{(j)})$ and estimate
the centers of the spectral measure and their respective masses as

(5.11) $\hat{\bf
c}_{j}=\frac{1}{\hat{N}_{n}^{(j)}}\sum_{i\in\hat{\mathcal{I}}_{n}^{(j)}}\frac{{\bf
X}_{i}}{\|{\bf X}_{i}\|}\quad\mbox{ and
}\quad\hat{\pi}_{j}=\frac{\hat{N}_{n}^{(j)}}{N_{n}},\quad j=1,\dots,p.$

The following result is an inmediate consequence of the main results of this
section.

######  Corollary 5.9.

Suppose $m=p$ and that the conditions of Proposition 5.1 and Theorem 5.6 hold.
Then, $\hat{\bf c}_{j}\stackrel{{\scriptstyle P}}{{\to}}{\bf c}_{j}$ and
$\hat{\pi}_{j}\stackrel{{\scriptstyle P}}{{\to}}w^{-1}\|{\bf
a}^{(j)}\|^{\alpha}$ for all $j=1,\dots,p$.

Note that in practice one can also normalize the estimates $\hat{\bf c}_{j}$
to ensure that they lie in the unit sphere for all $n$. Clearly the resulting
estimators remain consistent under the conditions of Corollary 5.9.

Even though the theoretical results of this section use the assumption
$\alpha>1$ in (5.8), we believe they should also hold when $\alpha\in(0,1]$.
The numerical experiments shown in the next section supports this assertion.

##  6\. Numerical illustrations

In all the examples considered below we compute weighted adjacency matrices
using the exponential kernel $d({\bf x},{\bf y})=\exp(-s\|{\bf x}-{\bf y}\|)$
with $s=1$ and select the number of clusters as suggested by the screeplots of
the fully connected weighted adjacency matrices ${W}$. It matched nicely the
correct number of clusters, when known. We consider sample sizes
$n=\\{1000,5000,25000,125000\\}$ and take a sample of extremes corresponding
to observations whose Euclidean norm is larger or equal to the following
vector of corresponding sample quantiles: $\beta=\\{0.9,0.96,0.984,0.9968\\}$,
respectively. These quantiles were chosen to lead to samples of extremes of
sizes $N_{n}=\\{100,200,400,800\\}$, correspondingly. For these extremes we
define $k_{n}$-nearest neighbors graphs with ${k_{n}=\lceil\frac{N_{n}}{C\log
N_{n}}\rceil+1}$, the corresponding values of the constants are in the vector
$C=\\{3,5,7,9\\}$ .

###  6.1. Linear factor model with and without noise

As a first example, consider $d-$dimensional vectors that follow the
$p-$dimensional linear factor model

(6.1) ${\bf X}=A{\bf Z}+\sigma\bm{\varepsilon}\,,$

where $A\in\mathbb{R}^{d\times p}$ is a matrix of factor loadings, ${\bf
Z}=(Z_{1},\ldots,Z_{p})^{\top}$ is a $p$-dimensional vector consisting of iid
standard Fréchet distributed components ($\alpha=1$), $\sigma\geq 0$ regulates
the signal to noise ratio and $\bm{\varepsilon}$ is a noise vector obtained by
multiplying a univariate independent standard Fréchet with an independent
$p$-dimensional random vector of iid standard normals, i.e.,

(6.2) $\bm{\varepsilon}={\bf N}\eta\,,$

where $\eta$ is standard Fréchet, ${\bf N}=(N_{1},\ldots,N_{p})^{\top}$ is a
$p-$random vector consisting of iid standard normals, and ${\bf Z},\eta$, and
${\bf N}$ are independent. Now using computations similar to those given in
Section 4, it can be shown that

(6.3) $\displaystyle\frac{\mathbb{P}(\|{\bf X}\|>x)}{\mathbb{P}(Z_{1}>x)}$
$\displaystyle\sim$ $\displaystyle\frac{\mathbb{P}(\sum_{i=1}^{p}\|{\bf
a}^{(i)}\|^{2}Z_{i}^{2}+\sigma^{2}\|{\bf
N}\|^{2}\eta^{2}>x^{2})}{\mathbb{P}(Z_{1}>x)}$ $\displaystyle\sim$
$\displaystyle\frac{\sum_{i=1}^{p}\mathbb{P}(\|{\bf
a}^{(i)}\|^{2}Z_{i}^{2}>x^{2})+\mathbb{P}(\sigma^{2}\|{\bf
N}\|^{2}\eta^{2}>x^{2})}{\mathbb{P}(Z>x)}$ $\displaystyle\to$
$\displaystyle\sum_{i=1}^{p}\|{\bf a}^{(i)}\|+\sigma\mathbb{E}\|{\bf
N}\|\,,~{}~{}~{}\mbox{as $x\to\infty$,}$

where the last line follows from an application of Breiman’s lemma, see
Breiman (1965). Taking this calculation one step further, we find that the
angular measure associated with the model is (6.1) (see (4.5)),

(6.4) $\Gamma(\cdot)=w^{-1}\left(\sum_{i=1}^{p}\|{\bf
a}^{(i)}\|\delta_{\frac{{\bf a}^{(i)}}{\|{\bf
a}^{(i)}\|}}(\cdot)+\sigma\mathbb{E}\|{\bf N}\|\delta_{\frac{{\bf N}}{\|{\bf
N}\|}}(\cdot)\right)\,,$

where $w=\sum_{i=1}^{p}\|{\bf a}^{(i)}\|+\sigma\mathbb{E}\|{\bf N}\|$. In
other words, $\Gamma$ has discrete mass points at $\frac{{\bf a}^{(i)}}{\|{\bf
a}^{(i)}\|}$ with probability $\|{\bf a}^{(i)}\|/w$, $i=1,\ldots,p$ and a
uniform distribution ${\bf N}/\|{\bf N}\|$ on $\mathbb{S}^{d-1}$ with
probability $\sigma\mathbb{E}\|{\bf N}\|/w$. This latter piece corresponds to
the noise component $\sigma\epsilon$. So the goal here is to identify the
discrete components of $\Gamma$ using our method when the model does not
strictly follow the LFM. Figure 2 shows pairwise scatter plots of the angular
components of extremes generated from a pure signal and a noisy LFM with
$\sigma>0$.

(a) Pure signal LFM (b) Noisy LFM

Figure 2. Pairwise scatterplots of the angular part of the extremes generated
from (6.1) with factor loading matrix (6.5), $n=125000$, $N_{n}=400$ and
$\sigma=\\{0,1\\}$. In both cases there are two clear clusters corresponding
to the signal. The red points in subfigure (b) denote extremes attributed to
the signal $A^{\top}\mathbf{Z}_{i}$.

We note that if $\sigma=0$, then model (6.1) is approximately equal to the
max-linear model
$X=(\vee_{j=1}^{k}a_{1j}Z_{j},\dots,\vee_{j=1}^{k}a_{pj}Z_{j})^{\top}$ and
will in fact have the same asymptotic spectral measure. Intuitively, this
model generates $p$ clusters of extremes since the noise term is only adding
uniform noise to the angular measure.

As part of a simulation study, we consider $\sigma=\\{0,1,3,5\\}$ and choose

(6.5) $A=\begin{pmatrix}0.1&0.9\\\ 0.2&0.8\\\ 0.3&0.7\\\ 0.4&0.6\\\
\end{pmatrix}.$

This model is similar to one of the max-linear models considered in the
simulations of Janßen and Wan (2020) where our factor loading matrix $A$ can
be viewed as a deterministic version of their random factor loadings. In the
simulations we took two clusters for the pure signal model where $\sigma=0$
and three clusters for the noisy model when $\sigma>0$ as these values are
suggested by the typical screeplots we observed; see Figure 3.

We compute the normalized columns of the $A$ matrix which correspond to the
location of the point masses of the spectral distribution (these are the ${\bf
c}_{k},\,k=1,2$ in (4.7)). After applying our method to a single realization
of size $n=125000$ with $N_{n}=400$, $k_{n}=\lceil\frac{400}{5\log
400}\rceil+1=15$, visualized in the pairwise scatter plots of Figure 2, we
obtained the estimates of the ${\bf c}_{k}$ represented in Figure 4. These
masses on the sphere are estimated by taking the mean of all members in each
of the identified clusters, seen in Figure 5, and then normalizing it to lie
on the unit sphere. The two panels in Figure 4 correspond to the cases of 2
clusters and no noise and two clusters with uniform noise. In the first plot,
the heat map does a good job in recreating the relatives size of the mass
locations. In the second panel, the first two columns of the matrix, also
reproduce the relative sizes of the columns (increasing in the first and
decreasing in the second) of the $A$ matrix. The third column corresponds to
the cluster of points that have not been assigned to either of the first two
clusters. As such they are essentially scattered uniformly around the unit
sphere but away from the locations of the point masses corresponding to the
first two columns. This is reflected in the third cluster having more negative
values as indicated by the softer (red colors) in the heat map.

Figure 3. Screeplots of fully connected kernel matrix of pure signal and noisy
linear factor models noise models.

(a) pure signal (b) uniform noise on the sphere

Figure 4. The heat maps show the estimated cluster centers based on the
cluster assignments displayed in Figure 5. The extremal sample corresponds to
four dimensional extremes generated from LFM given by (6.1) with loading
matrix (6.5) and $\sigma=0$ and $\sigma=1$ respectively. In both cases we took
$n=125000$, $N_{n}=400$ and $k_{n}=15.$

(a) Pure signal LFM (b) Noisy LFM

Figure 5. Cluster assignments output of spectral clustering applied to data
generated from (6.1) with $n=125000$, $N_{n}=400$ and $\sigma=\\{0,1\\}$. In
both cases spectral clustering used an extremal $15$-NN graph.

A small simulation study was conducted for this LFM model with and without
noise. Based on the screeplots, we used 2 clusters in the noiseless case and 3
clusters in the case with noise. The two normalized columns of the $A$ were
estimated and the boxplot of the estimation error measured in Frobenius norm
are displayed in Figure 6 for $\sigma=0$ and in Figure 7 for the case
$\sigma>0$. The succession of boxplots in each row correspond to an increasing
$N_{n}$, with the centers and width becoming smaller. Note that the scales on
the plots change across the row. The boxplots in blue correspond to our method
with difference choices of nearest neighbors as a function of $C$, and the
yellow boxplot is based on the spherical $k$-means approach considered in
Janßen and Wan (2020). In the $\sigma=0$ case, our method performs about the
same or slightly better than the spherical $k$-means method. However, as one
adds noise to the model, our method generally outperforms spherical $k$-means.
In models with larger noise, it can be more difficult to estimate the LFM
signal. So to compare performance across difference sample sizes and level of
noise, we can calibrate by calculating a notion of signal to noise ratio. In
this context we consider the part of the mass in the angular measure
associated to the signal in (6.1), which as a function of $\sigma$ is given by

$\mbox{SNR}(\sigma):=\frac{\sum_{i=1}^{p}\|{\bf
a}^{(i)}\|}{\sum_{i=1}^{p}\|{\bf a}^{(i)}\|+\sigma\mathbb{E}\|{\bf N}\|}\,.$

In the absence of any noise, i.e., $\sigma=0$, then SNR is 1 while as
$\sigma\to\infty$, SNR converges to 0. For the simulation example above for
which $d=4,p=2$, we have $\mathbb{E}\|{\bf
N}\|=\sqrt{2}\Gamma(5/2)/\Gamma(2)=1.880$. Hence
SNR$(\sigma)=2.065/(2.065+\sigma 1.880)$. In looking at the various plots in
Figure 7, it is instructive to compute the effective sample size given by
ESS$=\mbox{SNR}\times N_{n}$. This number essentially gives the expected
sample size of the number of observations, from the total $N_{n}$, that come
from the signal. With this index in mind, plots that have the same ESS values
(reported in the caption of Figure 7) generally show similar results since the
procedures are applied to the roughly the same number of extreme observations
attributed to the signal component in the model. We finally note that in this
simulation $\alpha=1$ which is not currently covered by our LFM theory but is
the setting proposed in the simulations of Janßen and Wan (2020). We carried
out simulations with $\alpha=0.5$ and $\alpha=2$ and obtained qualitatively
the same type of results as the ones reported here. The only noticeable
difference was that spherical $k$-means seems to work better with larger
$\alpha$ in the noisy model, but is much worse for small $\alpha$. In both
cases spectral clustering outperformed spherical $k$-means.

Figure 6. Estimation error measured in Frobenius norm when $\sigma=0$.

Figure 7. Estimation error of cluster centers measured in Frobenius norm. data
was generated from the noisy LFM (6.1). The sample sizes increases from left
to right as $n=\\{1000,5000,25000,125000\\}$ and
$N_{n}=\\{100,200,400,800\\}$, and from noise level increases from the top
down as $\sigma=\\{1,3,5\\}$. Across rows the ESS are: $\\{52,105,209,418\\}$,
$\\{27,54,107,214\\}$, $\\{18,36,72,144\\}$

###  6.2. Bivariate extremes from MA(3)

We consider the model discussed in the introduction and represented in Figure
1. More specifically, the model is
$Y_{t}=Z_{t}+.5Z_{t-1}-.6Z_{t-2}+1.5Z_{t-3}$, where $\\{Z_{t}\\}$ is an iid
symmetric stable random variables with index $\alpha=1.8$. We analyze the
extremal dependence structure of the bivariate vector ${\bf
X}_{t}=(Y_{t},Y_{t-1})^{\top}$ by looking for clusters in the extremes of
${\bf X}_{t}$. This model can be written in the form (4.1) since we can define
${\bf Z}_{t}=(Z_{t},Z_{t-1},Z_{t-2},Z_{t-3},Z_{t-4})$ and hence

${\bf X}_{t}=\begin{pmatrix}1&0.5&-0.6&1.5&0\\\
0&1&0.5&-0.6&1.5\end{pmatrix}{\bf Z}_{t}.$

Note that even though in this case the sample $\\{{\bf X}_{t}\\}$ is not
independent, the asymptotic distribution obtained in Theorem 4.1 still holds.
In particular, the angular distribution is supported in the points (4.7) i.e.

$\displaystyle{\bf c}_{1,\pm}=\pm(1,0),\quad{\bf
c}_{2,\pm}=\pm(.5,1)/\sqrt{1.25},\quad{\bf
c}_{3,\pm}=\pm(-0.6,0.5)/\sqrt{0.61},$ $\displaystyle\quad{\bf
c}_{4,\pm}=\pm(1.5,-0.6)/\sqrt{2.61}\quad\mbox{and}\quad{\bf
c}_{5,\pm}=\pm(0,1).$

Figure 8. Screeplot of kernel matrix and clustering performance of 2
dimensional MA$(3)$ extremes when $n=25000$ and $N_{n}=400$. Figure 9.
Estimation error of the matrix of atoms of the spectral measure of the
symmetric MA(3) model. The sample sample sizes were
$n=\\{1000,5000,25000,125000\\}$ giving $N_{n}=\\{100,200,400,800\\}$. We used
$k_{n}=\lceil\frac{N_{n}}{C\log N_{n}}\rceil+1$ nearest neighbor graphs.

Figure 8 illustrates the behavior of spectral clustering for this model when
$N_{n}=400$ and $k_{n}=\lceil\frac{400}{2\log 400}\rceil+1=35$. It is worth
noting that in Figure 1 we had a larger sample size of $100,000$ and stricter
quantile threshold of $0.998$ resulting in smaller number of observations
considered as extremes, but with an empirical distribution visibly closer to
the prescribed asymptotic discrete distribution. Therefore the simulation
scenario considered here is more difficult. Figure 9 illustrates the
convergence of the method. While the spectral $k$-means method of Janßen and
Wan (2020) performs slightly better than our spectral clustering for
$N_{n}\leq 200$, our proposed method appears better with much smaller
variability for a larger number of extremes. The choice $k_{n}$ of nearest
neighbors did not appreciably impact the performance of spectral clustering
across the different sample sizes.

###  6.3. Air pollution data

We revisit the data analyzed by Heffernan and Tawn (2004) and Janßen and Wan
(2020). It is available in the R package texmex and consists of daily
measurements of five air pollutants in the city of Leeds, UK. It was collected
between 1994 and 1998, and split into summer and winter months yielding a
total of $578$ and $532$ observations respectively. Following standard
practice in multivariate extremes data analysis we standardize the marginal
distribution of the data to focus on the extremal dependence. More
specifically, we transform the marginals of the original observations
$\mathbf{X}_{i}$ as in Janßen and Wan (2020) i.e., we let

$Y_{ij}:=1/\\{1-F_{nj}(X_{ij})\\},$

where $F_{nj}(x)=\frac{1}{n}\sum_{i=1}^{n}\mathbbm{1}(X_{i}\leq x)$ denotes
the $j$th marginal empirical cumulative distribution function,
$x\in\mathbb{R}$ and $j=1,\dots,d$. We then proceed to define the extremal
observations as the $10$% of the transformed observations
$\\{\mathbf{Y}_{i}\\}$ with largest Euclidean norm and analyze their angular
components with our algorithm. We analyze this data using spectral clustering
with the exponential kernel and $s=1$ as in the simulated data. The screeplots
in Figure 10 suggest that one should consider 5 clusters for this data.

Figure 10. Screeplots of the kernel matrix of the air pollution data extremes
obtained with the exponential kernel and bandwidth parameter $s=1$ .

Figure 11 shows the estimated cluster centers $\mathbf{c}_{j}$ for
$j=1,\dots,5$. We note that the “elbow plot” considered by Janßen and Wan
(2020) suggested the authors to use 4 or 5 clusters in their article. Our
results for 5 clusters is consistent with their analysis. Specifically, the
normalized cluster centers in the heat map of Figure 11 show that the extremes
of the five air pollutants act mostly independent. Looking a bit more closely,
both NO and NO2 share common strength in clusters 2 and 3, which is much
stronger in winter than in summer. PM10 also shares a common source (cluster
2) with NO and NO2, which is more pronounced in winter than summer. For the O3
and NO2 pollutants, we examined time lagged dependence by applying the
spectral clustering algorithm to the vector ${\bf
X}_{t}=(X_{t},X_{t-1},X_{t-2},X_{t-3})^{T}$, where $X_{t}$ represents either
the measured value of O3 or NO2 on day $t$. The resulting heat plots for the
cluster centers (4) are displayed in Figures 12 (O3) and 13 (NO2). The super
and sub diagonals reflect some extremal dependence at time lag 1 for O3 in
both summer and winter. This dependence mostly dissipates after one day. The
situation for NO3 is a bit more complex. One still discerns some extremal
dependence at a one day lag as indicated by the high-temperature in the heat
maps along the diagonal and subdiagonal. However, some clusters have similar
shading for its center of mass, e.g., clusters 1 and 4 for winter, which
suggests poor delineation between the clusters. In addition, there is a
stronger day effect in the summer than winter for NO2 and the dependence does
not necessarily die out after one day lag as in the O3 case.

(a) summer (b) winter

8

Figure 11. Five dimensional extremes from air pollution summer and winter
data. The heat maps show the estimated cluster centers using spectral
clustering with 5 clusters and 9-nearest neighbors.

(a) summer (b) winter

Figure 12. Four dimensional time series data constructed with lags 0-3 of O3
for summer and winter data respectively. The heat maps show the estimated
cluster centers using spectral clustering with 5 clusters and 9-nearest
neighbors.

(a) summer (b) winter

Figure 13. Four dimensional time series data constructed with lags 0-3 of NO2
for summer and winter data respectively. The heat maps show the estimated
cluster centers using spectral clustering with 5 clusters and 9-nearest
neighbors.

##  7\. Discussion

In this work we introduced a spectral clustering approach for learning the
angular measure of multivariate extremes. We proved that this approach leads
to consistent clustering for a natural linear factor model and showed the good
finite sample performance of our methods in numerical experiments. The
encouraging results suggest the method might be applied in more general
contexts. We are particularly interested in exploring two type of extensions.
First, high dimensional scenarios where the dimension of the extremes $d$
might be larger than the number of observed extremes $N_{n}$. This would
require introducing appropriate notions of sparsity and regularization.
Second, it seems natural to investigate generative models that lead to
continuous angular measures in the limit. This scenario implies one would need
to carefully introduce more general definitions of extremal clusters and
different analysis of the convergence of $k$-nearest neighbor graphs.

## Appendix

Before proving Lemmas 5.4 and 5.5 we will give a result regarding random
partitions of uniform random variables that we will leverage as the continuity
of $F_{T_{n}}$ implies that $F_{T_{n}}(T_{n})\sim\mbox{Unif}(0,1)$. We remind
the reader that in our proofs $c>0$ represents a finite and non-zero constant
whose value may change from line-to-line.

######  Lemma 7.1.

Let $U_{1},\dots,U_{N}\overset{iid}{\sim}\mbox{Unif}(0,1)$ and consider the
random partition of the unit interval
$I_{j,N}=[\frac{j-1}{m_{N}},\frac{j}{m_{N}})$, where
$m_{N}=\frac{N}{\tau\log(N)}$, $\tau>1$ and $j=1,\dots,m_{N}$. Then, with
probability at least $1-\frac{N^{1-\tau}}{\tau\log(N)}{(1+N^{-0.2\tau})}$

  1. (i)

Every $I_{j,N}$ contains at least one of the variables $U_{1},\dots,U_{N}$.

  2. (ii)

No $I_{j,N}$ contains more than $3\tau\log(N)$ of the variables
$U_{1},\dots,U_{N}$.

###### Proof.

Consider the event $E_{N,1}=\\{\mbox{Every }I_{j,N}\mbox{ contains at least
one of the variables }U_{1},\dots,U_{N}\\}$ and note that a union bound gives

$\displaystyle\mathbb{P}(E_{N,1})$ $\displaystyle\geq
1-\sum_{j=1}^{m_{N}}\mathbb{P}(U_{k}\notin I_{j,N},\,\forall k=1,\dots,N)$
$\displaystyle=1-m_{N}\left(1-\frac{1}{m_{N}}\right)^{N}$ $\displaystyle\geq
1-m_{N}e^{-N/m_{N}}$ (7.1) $\displaystyle=1-\frac{N^{1-\tau}}{\tau\log(N)}\,.$

Now consider the event $E_{N,2}=\\{\mbox{No }I_{j,N}\mbox{ contains more than
}3\tau\log(N)\mbox{ of the variables }U_{1},\dots,U_{N}\\}$. It follows again
from a union bound that yields

$\displaystyle\mathbb{P}(E_{N,2})$ $\displaystyle\geq
1-\sum_{j=1}^{m_{N}}\mathbb{P}(I_{j,N}\mbox{ has more than }3\tau\log(N)\mbox{
of the }U_{1},\dots,U_{N})$ $\displaystyle=1-m_{n}\mathbb{P}(S_{N}>3\tau\log
N),$

where $S_{N}\sim\mbox{Bin}(N,\frac{1}{m_{N}})$. Invoking Bernstein’s
inequality see that

(7.2) $\mathbb{P}(E_{N,2})\geq{1-m_{N}e^{-\frac{1}{2}\frac{(3\tau\log
N-\tau\log N)^{2}}{\tau\log N+2\tau\log
N/3}}=1-\frac{N^{-(1.2\tau-1)}}{\tau\log N}}.$

Combining (Proof.) and (7.2) shows that $(i)$ and $(ii)$ hold with the desired
probability. ∎

### Proof of Lemma 5.4

We prove the lemma for positive $Z$. The same type of arguments work in the
symmetric case and are therefore omitted. Note that we can write for $y>0$,

$\displaystyle f_{T_{n}}(t)=\frac{w_{1}^{2}}{c_{p}u_{n}}\cdot$
$\displaystyle\int_{-\infty}^{\infty}\cdots\int_{-\infty}^{\infty}\Bigl{[}z_{1}f_{Z}(z_{1})\,f_{Z}(z_{2})\cdots
f_{Z}(z_{p-1})f_{Z}\bigl{[}\bigl{(}tz_{1}w_{1}^{2}/u_{n}-(c_{2}z_{2}+\cdots+c_{p-1}z_{p-1})\bigr{)}/c_{p}\bigr{]}$
$\displaystyle\cdot\mathbf{1}\Bigl{(}z_{1}>u_{n}/w_{1}^{1/\alpha},\,(a_{11}z_{1}+a_{12}z_{2}+\cdots+a_{1p}\tilde{z}_{p})^{2}+(a_{21}z_{1}+a_{22}z_{2}\cdots+a_{2p}\tilde{z}_{p})^{2}>u_{n}^{2}\Bigr{)}\Bigr{]}$
(7.3) $\displaystyle\hskip 397.48499ptdz_{1}\cdots dz_{p-1}$
$\displaystyle\div\left[\mathbb{P}\Bigl{(}(a_{11}Z_{1}+\cdots+a_{1p}Z_{p})^{2}+(a_{21}Z_{1}+\cdots+a_{2p}Z_{p})^{2}>u_{n}^{2},\
\ {Z_{1}>u_{n}/w_{1}^{1/\alpha}}\Bigr{)}\right]$
$\displaystyle:=M_{n}(t)/D_{n}\,,$

where
$\tilde{z}_{p}=(tz_{1}w_{1}^{2}/u_{n}-(c_{2}z_{2}+\cdots+c_{p-1}z_{p-1})\bigr{)}/c_{p}$,
$c_{i}=a_{2i},\,i=1,\ldots,p$. We already know that

(7.4) $D_{n}\sim cu_{n}^{-\alpha},\ n\to\infty.$

Next, from (5.7), $\sup_{z}f_{Z}(z)=M<\infty$, we conclude by (5.8) that

$\displaystyle
M_{n}(t)\leq{\frac{Mw_{1}^{2}}{c_{p}u_{n}}}\int_{u_{n}/w_{1}^{1/\alpha}}^{\infty}z_{1}f_{Z}(z_{1})\,dz_{1}\sim
cu_{n}^{-\alpha},\quad\mbox{ as }n\to\infty\,.$

Hence there exists an ${G}\in(0,\infty)$ such that for all $n$ large enough,

(7.5) $f_{T_{n}}(t)\leq{G}\ \ \text{for all $t$.}$

This shows $(i)$. Let us now turn to claim $(ii)$ for concreteness consider
$0<t\leq 1$. Note that, for large $n$, the indicator in (7.3) is bounded from
below by the indicator of the set $E=\\{C^{-1}u_{n}<z_{1}<Cu_{n},\,|z_{i}|\leq
1,\,i=2,\ldots,p-1\\}$ for some large $C$. Then, on $E$, the argument of the
last function $f_{Z}$ in (7.3) is within a compact interval, so we obtain

$\displaystyle M_{n}(t)\geq$ $\displaystyle
c{\frac{Mw_{1}^{2}}{c_{p}u_{n}}}\int_{C^{-1}u_{n}}^{Cu_{n}}z_{1}f_{Z}(z_{1})\,dz_{1}\sim
cu_{n}^{-\alpha}\,,$

where the last relation follows from a direct application of (5.8). Along with
(7.4) this establishes $(ii)$.

Finally, note that

$\displaystyle M_{n}(t)\leq$
$\displaystyle{\frac{w_{1}^{2}}{c_{p}u_{n}}}\int_{u_{n}/{w_{1}^{1/\alpha}}}^{\infty}z_{1}f_{Z}(z_{1})\int_{\mathbb{R}}\cdots\int_{\mathbb{R}}f_{Z}(z_{2})\cdots
f_{Z}(z_{p-1})\Bigl{[}$
$\displaystyle~{}~{}f_{Z}\bigl{[}\bigl{(}tz_{1}{w_{1}^{2}}/u_{n}-(c_{2}z_{2}+\cdots+c_{p-1}z_{p-1})\bigr{)}/c_{p}\bigr{]}\Bigr{]}dz_{1}\cdots
dz_{p-1}$ $\displaystyle=$
$\displaystyle\frac{{w_{1}^{2}}u_{n}}{c_{p}}\int_{u_{n}/{w_{1}^{1/\alpha}}}^{\infty}z_{1}f_{Z}(u_{n}z_{1})\int_{\mathbb{R}}\cdots\int_{\mathbb{R}}f_{Z}(z_{2})\cdots
f_{Z}(z_{p-1})\Bigl{[}$
$\displaystyle~{}~{}f_{Z}\bigl{[}\bigl{(}tz_{1}{w_{1}^{2}}-(c_{2}z_{2}+\cdots+c_{p-1}z_{p-1})\bigr{)}/c_{p}\bigr{]}\Bigr{]}dz_{1}\cdots
dz_{p-1}\,.$

Using the upper bound in (5.8) it is easy to see that for some $c>0$ and
sufficiently large $t$,

(7.6) $\displaystyle\int_{\mathbb{R}}\cdots\int_{\mathbb{R}}f_{Z}(z_{2})\cdots
f_{Z}(z_{p-1})\Bigl{[}f_{Z}\bigl{[}\bigl{(}t-(c_{2}z_{2}+\cdots+c_{p-1}z_{p-1})\bigr{)}/c_{p}\bigr{]}\Bigr{]}dz_{2}\cdots
dz_{p-1}$ $\displaystyle\leq ct^{-(\alpha+1)}\,.$

Indeed, the integral is, up to a constant, equal to the density of a linear
combination of $Z_{1},\ldots,Z_{p-1}$. Therefore, for all $y$ large enough,
uniformly in $n$,

$\displaystyle M_{n}(t)\leq$ $\displaystyle
cu_{n}\int_{1/{w_{1}^{1/\alpha}}}^{\infty}z_{1}f_{Z}(u_{n}z_{1})(tz_{1})^{-(\alpha+1)}\,dz_{1}\leq
cu_{n}^{-\alpha}t^{-(\alpha+1)}\,,$

where once again we have used the upper bound in (5.8). Together with (7.4)
this shows the upper bound in $(iii)$. The lower bound in $(iii)$ can be
established in an identical way using the lower bound in (5.8). $\Box$

### Proof of Lemma 5.5

It follows from Lemma 7.1 that, outside of an event $\Omega_{n}^{(1)}$ with
$\mathbb{P}\bigl{(}\Omega_{n}^{(1)}\bigr{)}\to 0$, each one of the intervals
$I_{i,n}$ contains at least one of the points

$T_{ni}=\frac{a_{21}Z_{2,i}^{(*,1)}+\dots+a_{p1}Z_{p,i}^{(*,1)}}{w_{1}^{2}Z_{1,i}^{(*,1)}/u_{n}},\quad
i=1,\ldots,N_{n}^{(1)},$

and none of the intervals contains more than $3\tau\log N_{n}^{(1)}$ of these
points. Note that (7.5) implies that

$\frac{\partial}{\partial
t}F_{T_{n}}^{-1}(t)=\frac{1}{f_{T_{n}}(t)}\geq\frac{1}{{G}},\quad\forall
t\in\mathbb{R}$

and hence by the fundamental theorem of calculus

$F_{T_{n}}^{-1}\left(\frac{i}{m_{n}}\right)-F_{T_{n}}^{-1}\left(\frac{i-1}{m_{n}}\right)\geq\frac{1}{{G}m_{n}}$

This shows that the length of the intervals $I_{i,n}$ satisfies

(7.7) $|I_{i,n}|\geq l_{n}/{G},\ \ i=1,\ldots,m_{n},$

where $l_{n}=\frac{1}{m_{n}}$. Since the conditional law of
$(Z_{1}/u_{n},Z_{2},\ldots,Z_{p})$ given (5.6) converges weakly, as
$n\to\infty$, to the law of

$\bigl{(}W_{\alpha}/w_{1},Z_{2},\ldots,Z_{p}\bigr{)}$

as defined in Theorem 4.1, we see that

$F_{T_{n}}\Rightarrow G:=\text{the law of}\ \
\frac{c_{2}Z_{2}+\cdots+c_{p}Z_{p}}{W_{\alpha}}.$

It follows that the values $F_{T_{n}}(t_{0})$ converge, as $n\to\infty$, to a
finite limit. Therefore, there is $0<\delta<1$ such that
$F_{T_{n}}^{-1}\bigl{(}(i-1)/m_{n}\bigr{)}\geq t_{0}$ for all $n$ large enough
and all $i\geq(1-\delta)m_{n}$. We conclude by Lemma 5.4 $(iii)$ that for such
$n$ and $i$,

(7.8) $\displaystyle l_{n}=$ $\displaystyle
F_{T_{n}}\Bigl{(}F_{T_{n}}^{-1}\bigl{(}i/m_{n}\bigr{)}\Bigr{)}-F_{T_{n}}\Bigl{(}F_{T_{n}}^{-1}\bigl{(}(i-1)/m_{n}\bigr{)}\Bigr{)}$
$\displaystyle\in$
$\displaystyle\bigl{(}D^{-1},D\bigr{)}\int_{F_{T_{n}}^{-1}\bigl{(}(i-1)/m_{n}\bigr{)}}^{F_{T_{n}}^{-1}\bigl{(}i/m_{n}\bigr{)}}t^{-(\alpha+1)}\,dt.$

Furthermore,

(7.9)
$\int_{F_{T_{n}}^{-1}\bigl{(}(i-1)/m_{n}\bigr{)}}^{F_{T_{n}}^{-1}\bigl{(}i/m_{n}\bigr{)}}t^{-(\alpha+1)}\,dt\geq\Bigl{(}F_{T_{n}}^{-1}\bigl{(}i/m_{n}\bigr{)}\Bigr{)}^{-(\alpha+1)}\Bigl{(}F_{T_{n}}^{-1}\bigl{(}i/m_{n}\bigr{)}-F_{T_{n}}^{-1}\bigl{(}(i-1)/m_{n}\bigr{)}\Bigr{)},$

while leveraging again Lemma 5.4 $(iii)$ we see that

(7.10)
$\frac{m_{n}-i}{m_{n}}=\int_{F_{T_{n}}^{-1}\bigl{(}i/m_{n}\bigr{)}}^{\infty}f_{T_{n}}(t)\,dt\leq
D\int_{F_{T_{n}}^{-1}\bigl{(}i/m_{n}\bigr{)}}^{\infty}t^{-(\alpha+1)}\,dt=\frac{D}{\alpha}\Bigl{(}F_{T_{n}}^{-1}\bigl{(}i/m_{n}\bigr{)}\Bigr{)}^{-\alpha}.$

Combining (7.9) and (7.10), we conclude that

$\int_{F_{T_{n}}^{-1}\bigl{(}(i-1)/m_{n}\bigr{)}}^{F_{T_{n}}^{-1}\bigl{(}i/m_{n}\bigr{)}}t^{-(\alpha+1)}\,dt\geq
c\left(\frac{m_{n}-i}{m_{n}}\right)^{(\alpha+1)/\alpha}\Bigl{(}F_{T_{n}}^{-1}\bigl{(}i/m_{n}\bigr{)}-F_{T_{n}}^{-1}\bigl{(}(i-1)/m_{n}\bigr{)}\Bigr{)},$

and so by (7.8),

$l_{n}\geq c\left(\frac{m_{n}-i}{m_{n}}\right)^{(\alpha+1)/\alpha}|I_{i,n}|.$

Since an upper bound can be obtained in the same way, we conclude that for
some $D_{1}\geq 1$, for all $n$ large enough and all $i\geq(1-\delta)m_{n}$,

(7.11)
$D_{1}^{-1}{l_{n}}\left(\frac{m_{n}-i}{m_{n}}\right)^{-(\alpha+1)/\alpha}\leq|I_{i,n}|\leq
D_{1}{l_{n}}\left(\frac{m_{n}-i}{m_{n}}\right)^{-(\alpha+1)/\alpha}.$

Choose $K>2D_{1}^{2}$, and choose $K_{0}$ so that

(7.12)
$1+\frac{K}{K_{0}}<\left(\frac{K}{2D_{1}^{2}}\right)^{\alpha/(\alpha+1)}.$

Consider an interval $I_{i,n}$ with $(1-\delta)m_{n}\leq i\leq m_{n}-K_{0}$.
It follows from (7.11) and the choice of $i$ that any point in $I_{i,n}$ is
closer to any point in $I_{i+1,n}$ than to any point in an interval $I_{j,n}$
with $j<i-K$. To see this it suffices to show that

(7.13) $|I_{i,n}|+{|I_{i+1,n}|<\sum_{k=1}^{K}|I_{i-k,n}|}$

Using (7.11) and $i\leq m_{n}-K_{0}$, the left hand side of (7.13) can be
upper bounded by

$\displaystyle|I_{i,n}|+{|I_{i+1,n}|}$ $\displaystyle\leq
D_{1}{l_{n}}\left(\frac{K_{0}}{m_{n}}\right)^{-(\alpha+1)/\alpha}+D_{1}{l_{n}}\left(\frac{K_{0}+1}{m_{n}}\right)^{-(\alpha+1)/\alpha}$
(7.14)
$\displaystyle<2D_{1}{l_{n}}\left(\frac{K_{0}}{m_{n}}\right)^{-(\alpha+1)/\alpha}.$

For $i\geq(1-\delta)m_{n}$ and using (7.11) once again, the right hand side of
(7.13) can be lower bounded by

$\displaystyle\sum_{k=1}^{K}{|I_{i-k,n}|}$ $\displaystyle\geq
D_{1}^{-1}{l_{n}}\sum_{k=1}^{K}\left(\frac{m_{n}-(i-k)}{m_{n}}\right)^{-(\alpha+1)/\alpha}$
$\displaystyle\geq
D_{1}^{-1}{l_{n}}K\left(\frac{m_{n}-i+K}{m_{n}}\right)^{-(\alpha+1)/\alpha}$
$\displaystyle\geq D_{1}^{-1}{l_{n}}K\left(\frac{\delta
m_{n}+K}{m_{n}}\right)^{-(\alpha+1)/\alpha}$ (7.15)
$\displaystyle>D_{1}^{-1}{l_{n}}K\left(1+\frac{K}{K_{0}}\right)^{-(\alpha+1)/\alpha}$

It follows from (Proof of Lemma 5.5) and (Proof of Lemma 5.5) that a suffcient
condition for establishing (7.13) is

$2D_{1}\left(\frac{K_{0}}{m_{n}}\right)^{-(\alpha+1)/\alpha}\leq
D_{1}^{-1}K\left(1+\frac{K}{K_{0}}\right)^{-(\alpha+1)/\alpha}.$

The last condition implies (7.12) for large $n$. We conclude that, on the
event $\Omega_{n}^{(1)}$, in a $k_{n}$-NN graph with

(7.16) $k_{n}>3(K+1)\tau\log N_{n}^{(1)},$

then all points $\bigl{(}V^{(j)},\,j=1,\ldots,N_{n}^{(1)}\bigr{)}$ within
$I_{i,n}$ in the range $(1-\delta)m_{n}\leq i\leq m_{n}-K_{0}$ will be
connected both to each other and to such a point in each $I_{i-1,n}$ and
$I_{i+1,n}$. The next observation to make is that, as long as $\delta$ is
small enough, the sequence $\bigl{(}F_{T_{n}}^{-1}(1-\delta)\bigr{)}$ is
bounded from above. Therefore, by Lemma 5.4 $(ii)$, uniformly in large enough
$n$, the density $f_{T_{n}}$ is bounded from below by, say, $a>0$ on the
interval $\bigl{(}0,F_{T_{n}}^{-1}(1-\delta)\bigr{)}$. Therefore, for all
large enough $n$,

(7.17) $|I_{i,n}|\,\leq\,l_{n}/a,\ \ 1\leq i\leq(1-\delta)m_{n}.$

To see this, note that

$\frac{\partial}{\partial
t}F_{T_{n}}^{-1}(t)=\frac{1}{f_{T_{n}}(t)}\leq\frac{1}{a},\quad\forall
t\in\mathbb{R}$

and hence by the fundamental theorem of calculus

$|I_{i,n}|=F_{T_{n}}^{-1}\left(\frac{i}{m_{n}}\right)-F_{T_{n}}^{-1}\left(\frac{i-1}{m_{n}}\right)\leq\frac{1}{am_{n}}=\frac{l_{n}}{a}$

It follows from (7.7) and (7.17) that if $K>\frac{2{G}}{a}$ then any point in
$I_{i,n}$ is closer to any point in $I_{i-1,n}$ and in $I_{i+1,n}$ than to any
point in an interval $I_{j,n}$ with $j<i-K$ or with $j>i+K$. Therefore, on the
event $\Omega_{n}^{(1)}$, in a $k_{n}$-NN graph satisfying (7.16), all points
$\bigl{(}T_{nj},\,j=1,\ldots,N_{n}^{(1)}\bigr{)}$ within $I_{i,n}$ in the
range $1\leq i\leq(1-\delta)m_{n}$ will be connected both to each other and to
such a point in each $I_{i-1,n}$ and $I_{i+1,n}$. Indeed, to show this it
suffices to show again that (7.13) holds true in the range $1\leq
i\leq(1-\delta)m_{n}$. It is easy to see that (7.7), (7.17) and
$K>\frac{2{G}}{a}$ entail

$|I_{i,n}|+{|I_{i+1,n}|}\leq\frac{2l_{n}}{a}<\frac{Kl_{n}}{{G}}\leq\sum_{k=1}^{K}{|I_{i-k,n}|.}$

Finally, it is obvious that if $K>K_{0}$, then on the same event
$\Omega_{n}^{(1)}$, in a $k_{n}$-NN graph satisfying (7.16), all points
$\bigl{(}T_{nj},\,j=1,\ldots,N_{n}^{(1)}\bigr{)}$ within $I_{i,n}$ in the
range $m_{n}-K_{0}<i\leq m_{n}$ will be connected both to each other and to a
such a point in each $I_{i-1,n}$ and $I_{i+1,n}$.

Summarizing the above discussion we conclude that on the event
$\Omega_{n}^{(1)}$, in a $k_{n}$-NN graph satisfying (7.16) with $K$ large
enough, all points $\bigl{(}T_{nj},\,j=1,\ldots,N_{n}^{(1)}\bigr{)}$ within
$I_{i,n}$ in the entire range $1\leq i\leq m_{n}$ will be connected both to
each other and to a such a point in each $I_{i-1,n}$ and $I_{i+1,n}$. In
particular, the $k_{n}$-NN graph will be connected.

We now translate this discussion to the random vectors ${\bf
M}^{(i)},\,i=1,\ldots,N_{n}^{(1)}$. We define intervals along vector ${\bf b}$
by

$J_{i,n}=I_{i,n}{\bf b},\ =1,\ldots,N_{n}^{(1)}\,.$

Then, outside of the event $\Omega_{n}^{(1)}$, each one of these intervals
contains at least one of the points $\bigl{(}{\bf
M}^{(i)},\,i=1,\ldots,N_{n}^{(1)}\bigr{)}$ and none of the intervals contains
more than $3\tau\log N_{n}^{(1)}$ of these points. By (7.7) the lengths of
these intervals satisfy for some ${G}_{1}>0$,

$|J_{i,n}|\geq l_{n}/{G}_{1},\ \ i=1,\ldots,m_{n}\,.$

We finally note that by (4.18), with probability tending to one $N_{n}\sim
Cnu_{n}^{-\alpha}$, and therefore ${G}>0$ and $n$ large enough ensure that
(7.16) holds provided $k_{n}>{G}\log n$. This concludes the proof. $\Box$

### Proof of Theorem 5.6

Lemma 5.5 gives us the connectivity of the extremal $k_{n}$-NN graph for
$k_{n}$ satisfying (7.16) with $K$ large enough. The next step is to
understand by how much the points $\bigl{(}{\bf
M}^{(i)},\,i=1,\ldots,N_{n}^{(1)}\bigr{)}$ are shifted by adding to them
$\bigl{(}{\bf D}^{(i)},\,i=1,\ldots,N_{n}^{(1)}\bigr{)}$ in (5.4). Denote
$\Omega_{n}^{(2)}={B}_{n}^{c}$ as defined in Lemma 4.4. Then
$\mathbb{P}\bigl{(}\Omega_{n}^{(2)}\bigr{)}\to 0$ as $n\to\infty$ and it is
elementary to check that outside of $\Omega_{n}^{(2)}$ we have $\|{\bf
D}^{(i)}\|\leq ch_{n}^{2}/u_{n}$ for all $i=1,\ldots,N_{n}^{(1)}$. Recall that
by the choice of $h_{n}$ we have

$h_{n}^{2}/u_{n}=o(l_{n})\ \ \text{as $n\to\infty$.}$

If we define new sets by

$\tilde{J}_{i,n}=\bigl{\\{}{\bf M}^{(j)}+{\bf D}^{(j)}:\,{\bf M}^{(j)}\in
J_{i,n}\bigr{\\}},\ i=1,\ldots,m_{n},$

then it follows immediately that for large $n$, outside of the event
$\Omega_{n}^{(2)}$, the new sets have the property described by Lemma 5.5,
perhaps with a larger $K_{0}$. We already know that this means that for large
$n$, outside of $\Omega_{n}^{(1)}\cup\Omega_{n}^{(2)}$, the extremal
$k_{n}$-NN graph with $k_{n}$ satisfying (7.16) with $K$ large enough, is
connected. $\Box$

## References

  * Basrak et al. [2002] Bojan Basrak, Richard A Davis, and Thomas Mikosch.  A characterization of multivariate regular variation.  _The Annals of Applied Probability_ , 12(3):908–920, 2002.
  * Belkin and Niyogi [2003] Mikhail Belkin and Partha Niyogi.  Laplacian eigenmaps for dimensionality reduction and data representation.  _Neural computation_ , 15(6):1373–1396, 2003\.
  * Breiman [1965] Leonard Breiman.  On some limit theorems similar to the arc-sin law.  _Theory of Probability & Its Applications_, 10(2):323–331, 1965.
  * Chautru [2015] Emilie Chautru.  Dimension reduction in multivariate extreme value analysis.  _Electronic Journal of Statistics_ , 9(1):383–418, 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.  _arXiv preprint arXiv:2104.03966_ , 2021.
  * Cooley and Thibaud [2019] Daniel Cooley and Emeric Thibaud.  Decompositions of dependence for high-dimensional extremes.  _Biometrika_ , 106(3):587–604, 2019.
  * Davis and Resnick [1989] Richard A Davis and Sidney I Resnick.  Basic properties and prediction of max-arma processes.  _Advances in Applied Probability_ , 21(4):781–803, 1989.
  * Davison and Huser [2015] Anthony C Davison and Raphaël Huser.  Statistics of extremes.  _Annual Review of Statistics and its Application_ , 2:203–235, 2015.
  * Dhillon and Modha [2001] Inderjit S Dhillon and Dharmendra S Modha.  Concept decompositions for large sparse text data using clustering.  _Machine Learning_ , 42(1):143–175, 2001.
  * Drees and Sabourin [2021] Holger Drees and Anne Sabourin.  Principal component analysis for multivariate extremes.  _Electronic Journal of Statistics_ , 15(1):908–943, 2021.
  * Embrechts and Goldie [1980] Paul Embrechts and Charles M Goldie.  On closure and factorization properties of subexponential and related distributions.  _Journal of the Australian Mathematical Society_ , 29(2):243–256, 1980.
  * Engelke and Ivanovs [2021] Sebastian Engelke and Jevgenijs Ivanovs.  Sparse structures for multivariate extremes.  _Annual Review of Statistics and Its Application_ , 8:241–270, 2021.
  * Gissibl and Klüppelberg [2018] Nadine Gissibl and Claudia Klüppelberg.  Max-linear models on directed acyclic graphs.  _Bernoulli_ , 24(4A):2693–2720, 2018.
  * Goix et al. [2015] Nicolas Goix, Anne Sabourin, and Stéphan Clémen.  Learning the dependence structure of rare events: a non-asymptotic study.  In _Conference on Learning Theory_ , pages 843–860. PMLR, 2015.
  * 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.
  * Hall et al. [2002] Peter Hall, Liang Peng, and Qiwei Yao.  Moving-maximum models for extrema of time series.  _Journal of statistical planning and inference_ , 103(1-2):51–63, 2002.
  * Heffernan and Tawn [2004] Janet E Heffernan and Jonathan A Tawn.  A conditional approach for multivariate extreme values (with discussion).  _Journal of the Royal Statistical Society: Series B_ , 66(3):497–546, 2004.
  * Hendrickson and Leland [1995] Bruce Hendrickson and Robert Leland.  An improved spectral graph partitioning algorithm for mapping parallel computations.  _SIAM Journal on Scientific Computing_ , 16(2):452–469, 1995.
  * Janßen and Wan [2020] Anja Janßen and Phyllis Wan.  $k$-means clustering of extremes.  _Electronic Journal of Statistics_ , 14(1):1211–1233, 2020.
  * Klüppelberg and Lauritzen [2019] Claudia Klüppelberg and Steffen Lauritzen.  Bayesian networks for max-linear models.  In _Network Science_ , pages 79–97. Springer, 2019.
  * Lei and Rinaldo [2015] Jing Lei and Alessandro Rinaldo.  Consistency of spectral clustering in stochastic block models.  _The Annals of Statistics_ , 43(1):215–237, 2015\.
  * Meyer and Wintenberger [2019] Nicolas Meyer and Olivier Wintenberger.  Sparse regular variation.  _arXiv preprint arXiv:1907.00686_ , 2019.
  * Ng et al. [2002] Andrew Y Ng, Michael I Jordan, and Yair Weiss.  On spectral clustering: Analysis and an algorithm.  In _Advances in Neural Information Processing Systems_ , pages 849–856, 2002.
  * Resnick [2007] Sidney I Resnick.  _Heavy-tail phenomena: probabilistic and statistical modeling_.  Springer Science & Business Media, 2007.
  * Resnick [2018] Sidney I Resnick.  _Extreme values, regular variation and point processes_.  Springer, New York, 2018.
  * Rohe et al. [2011] Karl Rohe, Sourav Chatterjee, and Bin Yu.  Spectral clustering and the high-dimensional stochastic blockmodel.  _The Annals of Statistics_ , 39(4):1878–1915, 2011.
  * Shi and Malik [2000] Jianbo Shi and Jitendra Malik.  Normalized cuts and image segmentation.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 22(8):888–905, 2000.
  * Simpson et al. [2020] Emma S Simpson, Jennifer L Wadsworth, and Jonathan A Tawn.  Determining the dependence structure of multivariate extremes.  _Biometrika_ , 107(3):513–532, 2020.
  * Van Driessche and Roose [1995] Rafael Van Driessche and Dirk Roose.  An improved spectral bisection algorithm and its application to dynamic load balancing.  _Parallel Computing_ , 21(1):29–48, 1995.
  * Von Luxburg [2007] Ulrike Von Luxburg.  A tutorial on spectral clustering.  _Statistics and computing_ , 17(4):395–416, 2007\.
  * Zhou and Amini [2019] Zhixin Zhou and Arash A Amini.  Analysis of spectral clustering algorithms for community detection: the general bipartite setting.  _The Journal of Machine Learning Research_ , 20(1):1774–1820, 2019.