A Multivariate Extreme Value Theory Approach to Anomaly Clustering and Visualization

TL;DR

This paper introduces a novel multivariate extreme value theory-based mixture model for anomaly clustering and visualization, enabling better identification and interpretation of extreme events in complex systems.

Contribution

It proposes a new mixture model that treats anomaly types as latent variables, facilitating clustering and visualization of extreme observations.

Findings

Effective clustering of extreme anomalies demonstrated on simulated data.

Successful application to aeronautics data shows practical utility.

Posterior probabilities provide meaningful similarity measures between anomalies.

Abstract

In a wide variety of situations, anomalies in the behaviour of a complex system, whose health is monitored through the observation of a random vector X = (X1,. .. , X d) valued in R d , correspond to the simultaneous occurrence of extreme values for certain subgroups $\alpha$ $\subset$ {1,. .. , d} of variables Xj. Under the heavy-tail assumption, which is precisely appropriate for modeling these phenomena, statistical methods relying on multivariate extreme value theory have been developed in the past few years for identifying such events/subgroups. This paper exploits this approach much further by means of a novel mixture model that permits to describe the distribution of extremal observations and where the anomaly type $\alpha$ is viewed as a latent variable. One may then take advantage of the model by assigning to any extreme point a posterior probability for each anomaly type $\alpha$, defining implicitly a similarity measure between anomalies. It is explained at length how the latter permits to cluster extreme observations and obtain an informative planar representation of anomalies using standard graph-mining tools. The relevance and usefulness of the clustering and 2-d visual display thus designed is illustrated on simulated datasets and on real observations as well, in the aeronautics application domain.

Tables

Table 1. Table 1: Average error on the model parameters, n 0 = 1 ​ e ​ 3 subscript 𝑛 0 1 𝑒 3 n_{0}=1e3 extreme points

	${\lambda }_k=1.$	${\lambda }_k=0.75$	${\lambda }_k=0.5$	${\lambda }_k=0.25$	${\lambda }_k=0.1$
$err\left(\hat{\rho }\right)$	1.39e-5	1.37e-5	1.57e-5	1.22e-5	2.11e-5
$err\left(\hat{\nu }\right)$	5.53	5.81	6.28	6.41	9.06
$err\left(\hat{\lambda }\right)$	2.65e-2	2.04e-2	1.19e-2	5.97e-3	3.66e-3

Table 2. Table 2: Average error on the model parameters, n 0 = 2 ​ e ​ 3 subscript 𝑛 0 2 𝑒 3 n_{0}=2e3 extreme points

	${\lambda }_k=1.$	${\lambda }_k=0.75$	${\lambda }_k=0.5$	${\lambda }_k=0.25$	${\lambda }_k=0.1$
$err\left(\hat{\rho }\right)$	9.98e-6	1.12e-5	1.06e-5	1.62e-5	1.64e-5
$err\left(\hat{\nu }\right)$	3.23	4.13	4.08	4.29	5.05
$err\left(\hat{\lambda }\right)$	1.62e-2	1.2e-2	8.11e-3	4.28e-3	3.11e-3

Table 3. Table 3: Average number of labeling errors

	${\lambda }_k=1.$	${\lambda }_k=0.75$	${\lambda }_k=0.5$	${\lambda }_k=0.25$	${\lambda }_k=0.1$
$n_0=1e3$	0.	0.	0.	0.6	264.4
$n_0=2e3$	0.	0.	0.4	1.8	537.8

Table 4. Table 4: Purity score - Comparisons with standard approaches for different extreme sample sizes.

	$n_0=500$	$n_0=400$	$n_0=300$	$n_0=200$	$n_0=100$
Dirichlet mixture	0.8	0.82	0.82	0.84	0.85
Kmeans	0.72	0.73	0.75	0.78	0.8
Spectral clustering	0.78	0.77	0.82	0.81	0.8

Equations

$$n\mathbb{P}\left(n^{-1}\mathbf{V}\in A\right)\to \mu (A)$$

$$\mu \circ T^{-1}(\mathrm{d}r,\mathrm{d}\mathbf{w})=r^{-2}\mathrm{d}r\Phi (\mathrm{d}\mathbf{w}).$$

$$\mathbb{P}\left(\mathbf{W}\in A,R>tr|R>t\right)$$

$${\int }_{{\mathcal{S}}_d}{w}_i\Phi (\mathrm{d}\mathbf{w})=1,\text{ for }i=1,\dots ,d.$$

$$\Phi ({\mathcal{S}}_d)={\int }_{{\mathcal{S}}_d}\Phi (\mathrm{d}\mathbf{w})={\int }_{{\mathcal{S}}_d}(w_1+\dots +w_d)\Phi (\mathrm{d}\mathbf{w})=d.$$

$$\mathcal{C}_{\alpha}=\big{\{}\mathbf{v}\geq 0:~{}\|\mathbf{v}\|_{\infty}\geq 1,~{}v_{i}>0\text{ for }i\in\alpha,\;v_{i}=0~{}\text{ for }i\notin\alpha\big{\}}.$$

$$\mathcal{R}_{\alpha}^{\epsilon}~{}=~{}\big{\{}\mathbf{v}\geq 0,~{}\|\mathbf{v}\|_{\infty}\geq 1,~{}v_{i}>\epsilon~{}\text{ for }i\in\alpha,\;v_{i}\leq\epsilon~{}\text{ for }i\notin\alpha\big{\}},$$

$$\mu ({\mathcal{C}}_{\alpha })=\underset{ϵ\to 0}{\lim }\mu ({\mathcal{R}}_{\alpha }^{ϵ})$$

$$\mu ({\mathcal{R}}_{\alpha }^{ϵ})=\underset{t\to \infty }{\lim }t\mathbb{P}\left(\right)(\Vert \mathbf{V}{\Vert }_{\infty }>t,\forall j\in \alpha :V_j>tϵ,\forall j\notin \alpha :V_j<tϵ).$$

$$d^{-1}\Phi (\cdot )=\sum _{k=1}^K{\pi }_k{\Phi }_{{\alpha }_k}(\cdot )+\sum _{j\le d_1}{\pi }_{K+j}{\delta }_{{\mathbf{e}}_j}(\cdot ),$$

$$d^{-1}=\sum _{k=1}^K{\int }_{{\mathcal{S}}_{{\alpha }_k}}{w}_i{\pi }_k{\Phi }_{{\alpha }_k}(\mathrm{d}\mathbf{w})+\sum _{j\le d_1}{\int }_{{\mathcal{S}}_{\{j\}}}{w}_i{\pi }_{K+j}{\delta }_{{\mathbf{e}}_j}(\mathrm{d}\mathbf{w})={\pi }_{K+i}.$$

$$\Phi (\cdot )=d\sum _{k=1}^K{\pi }_k{\Phi }_{{\alpha }_k}(\cdot )+\sum _{j\le d_1}{\delta }_{{\mathbf{e}}_j}(\cdot ),$$

$$\sum _{k=1}^K{\pi }_k=1-d_1/d.$$

$${\phi }_{\alpha }(\mathbf{w}|{\mathbf{m}}_{\alpha },{\nu }_{\alpha })=\frac{\Gamma ({\nu }_{\alpha })}{{\prod }_{i\in \alpha }\Gamma ({\nu }_{\alpha }{m}_{\alpha ,i})}\prod _{i\in \alpha }{w}_i^{{\nu }_{\alpha }{m}_{\alpha ,i}-1}.$$

$$\frac{1}{d}$$

$$\mathbf{V}={\mathbf{V}}_k+{\mathbf{\epsilon }}_k=R_k{\mathbf{W}}_k+{\mathbf{\epsilon }}_k,$$

$$\{\begin{array}{rl}& \mathbb{P}\left(R_k>r\right)=r_0{r}^{-1},r>r_0,\\ & {\mathbf{W}}_k\sim {\Phi }_k,\\ & {\epsilon }_j\sim 1+\mathcal{E}xp({\lambda }_k),j\in \{1,\dots ,d\}\setminus {\alpha }_k,\end{array}$$

$${\gamma }_{i,k}^{(t+1)}=\mathbb{P}\left(\right)(Z_{i,k}=1|{\mathbf{V}}_i,{\mathbf{\theta }}_t)$$

$$Q(\mathbf{\theta },{\gamma }^{(t)})=\sum _i\sum _k{\gamma }_{i,k}^{(t+1)}\log p({\mathbf{V}}_i|Z_{i,k}=1,\mathbf{\theta }).$$

$$\begin{gathered} p(\mathbf{v}|\mathbf{\theta })=r_0\sum _{k=1}^K{\pi }_k{r}_k^{-|{\alpha }_k|-1}{\phi }_k({\mathbf{w}}_k|{\mathbf{m}}_k,{\nu }_k)\prod _{j\in {\alpha }_k^c}{f}_{\epsilon }(v_j|{\lambda }_k) \\ +\frac{r_0}{d}\sum _{k=K+1}^{K+d_1}{r}_k^{-2}\prod _{j\in \{1,\dots ,d\}\setminus k}{f}_{\epsilon }(v_j|{\lambda }_k) \end{gathered}$$

$${\nu }_k>0(1\le k\le K),{\lambda }_k>0(1\le k\le K+d_1),$$

$${\pi }_k:=\sum _{j\in {\alpha }_k}{\rho }_{k,j}\text{ and }{m}_{k,j}:=\frac{{\rho }_{k,j}}{{\pi }_k},\forall j\in {\alpha }_k.$$

$$\sum _{\{k:j\in {\alpha }_k\}}{\rho }_{k,j}$$

$$p(\mathbf{v}|z_k=1,\mathbf{\theta })=\{\begin{array}{l}{r}_k^{-|{\alpha }_k|-1}{\phi }_k({\mathbf{w}}_{\mathbf{k}}|{\rho }_{\mathbf{k}},{\nu }_k){\prod }_{j\in {\alpha }_k^c}{f}_{\epsilon }(v_j|{\lambda }_k),(k\le K)\\ v_k^{-2}{\prod }_{j\in \{1,\dots ,d\}\setminus k}{f}_{\epsilon }(v_j|{\lambda }_k),(K<k\le K+d_1).\end{array}$$

$${\pi }_k=\{\begin{array}{ll}{\sum }_{\ell \in {\alpha }_k}{\rho }_{k,l}& \text{ for }1\le k\le K,\\ d^{-1}& \text{ for }K<k\le K+d_1,\end{array}$$

$${\gamma }_{i,k}=p(Z_{i,k}=1|{\mathbf{V}}_i,\mathbf{\theta })=\frac{{\pi }_{k}p({\mathbf{V}}_i|z_{i,k}=1,\mathbf{\theta })}{{\sum }_{\begin{array}{c}1\le \ell \le K+d_1\\ \ell \ne k\end{array}}{\pi }_{\ell }p({\mathbf{V}}_i|z_{i,\ell }=1,\mathbf{\theta })},$$

$$Q(\bm{\theta},\bm{\gamma})=\sum_{i=1}^{n}\Big{[}\sum_{k=1}^{K}\gamma_{i,k}\big{[}\log\pi_{k}+\log\varphi_{k}(\mathbf{W}_{i,k}|\bm{\rho}_{k},\nu_{k})+\sum_{l\in\alpha_{k}^{c}}\log f_{\varepsilon}(V_{i,l}|\lambda_{k})\big{]}\\ +\sum_{k=K+1}^{K+d_{1}}\gamma_{i,k}\big{[}\sum_{\ell\neq k}\log f_{\varepsilon}(V_{i,l}|\lambda_{k})\big{]}\Big{]}+C=Q_{1}(\bm{\rho},\bm{\nu})+Q_{2}(\bm{\lambda})+C,$$

$$Q_1(\mathbf{\rho },\mathbf{\nu })$$

$$Q_2(\mathbf{\lambda })$$

$$\underset{\begin{array}{c}\mathbf{\theta }\text{ s.t.}\\ \text{eqref}eq:positiveConstraints,\text{eqref}eq:rho-constraint\end{array}}{\max }Q(\mathbf{\theta },\mathbf{\gamma })=\underset{\begin{array}{c}\mathbf{\rho },\mathbf{\nu }\text{ s.t.}\\ \text{eqref}eq:rho-constraint,{\nu }_k>0,k\le K\end{array}}{\max }{Q}_1(\mathbf{\rho },\mathbf{\nu })+\underset{\begin{array}{c}\mathbf{\lambda }\text{ s.t.}\\ {\lambda }_k>0,1\le k\le K+d_1\end{array}}{\max }{Q}_2(\mathbf{\lambda })$$

$${\lambda }_k^{\ast }=\frac{|{\alpha }_k^c|{\sum }_{i=1}^n{\gamma }_{i,k}}{{\sum }_{i=1}^n{\gamma }_{i,k}{\sum }_{l\in {\alpha }_k^c}(V_{i,\ell }-1)},k\le K+d_1.$$

Full text

Keywords— Anomaly detection, clustering, graph-mining, latent variable analysis, mixture modelling, multivariate extreme value theory, visualization

1 Introduction

Motivated by a wide variety of applications ranging from fraud detection to aviation safety management through the health monitoring of complex networks, data center infrastructure management or food risk analysis, unsupervised anomaly detection is now the subject of much attention in the data science literature, see e.g. Gorinevsky et al., (2012); T. Fawcett, (1997); Viswanathan et al., (2012). In frequently encountered practical situations and from the viewpoint embraced in this paper, anomalies coincide with rare measurements that are extremes, i.e. located far from central statistics such as the sample mean. In the $1$-d setting, numerous statistical techniques for anomaly detection are based on a parametric representation of the tail of the observed univariate probability distribution, relying on extreme value theory (EVT), see e.g. Clifton et al., (2011); Lee and Roberts, (2008); Roberts, (2000); Tressou, (2008) among others. In (even moderately) large dimensional situations, the modelling task becomes much harder. Many nonparametric heuristics for supervised classification have been adapted, substituting rarity for labeling, see e.g. Schölkopf et al., (2001), Steinwart et al., (2005) or Liu et al., (2008). In the unsupervised setting, several extensions of the basic linear Principal Component Analysis for dimensionality reduction and visualization techniques have been proposed in the statistics and data-mining literature, accounting for non linearities or increasing robustness for instance, cf Gorban et al., (2008) and Kriegel et al., (2008). These approaches intend to describe parsimoniously the ‘center’ of a massive data distribution, see e.g. Naik, (2017) and the references therein. Similarly, for clustering purposes, several multivariate heavy-tailed distributions have been proposed that are robust to the presence of outliers, see e.g. Forbes and Wraith, (2014), Punzo and Tortora, (2018). However the issue of clustering extremes or outliers is only recently receiving attention, at the instigation of industrial applications such as those mentioned above and because of the increasing availability of extreme observations in databases: generally out-of-sample in the past, extreme values are becoming observable in the Big Data era. It is the goal of the present article to propose a novel mixture model-based approach for clustering extremes in the multivariate setup, i.e. when the observed random vector $\mathbf{X}=(X_{1},\;\ldots,\;X_{d})$ takes its values in the positive orthant of the space $\mathbb{R}^{d}$ with $d>1$ equipped with the sum-norm $\|{(x_{1},\;\ldots,\;x_{d})}\|=\sum_{1\leq j\leq d}|x_{j}|$: ’extremes’ coinciding then with values $x$ such that $\mathbb{P}\left(\|{\mathbf{X}}\|>\|{\mathbf{x}}\|\right)$ is ’extremely small’. Precisely, it relies on a dimensionality reduction technique of the tail distribution recently introduced in Goix et al., (2017) and Goix et al., (2016), and referred to as the DAMEX algorithm. Based on multivariate extreme value theory (MEV theory), the latter method may provide a hopefully sparse representation of the support of the angular measure related to the supposedly heavy-tailed distribution of the random vector $\mathbf{X}$. As the angular measure asymptotically describes the dependence structure of the variables $X_{j}$ in the extremal domain (and, roughly speaking, permits to assign limit probabilities to directions $\mathbf{x}/\|{\mathbf{x}}\|$ in the unit sphere along which extreme observations may occur), this statistical procedure identifies the groups $\alpha\subset\{1,\;\ldots,\;d\}$ of feature indices such that the collection of variables $\{X_{j}:\;j\in\alpha\}$ may be simultaneously very large, while the others, the $X_{j}$’s for $j\notin\alpha$, remain small. Groups of this type are in 1-to-1 correspondence with the faces $\Omega_{\alpha}=\{\mathbf{x}\in\mathbb{R}^{d}:\;\|{\mathbf{x}}\|=1,\;x_{j}=0\text{ if }j\notin\alpha\text{ and }x_{j}>0\text{ if }j\in\alpha\}$ of the unit sphere composing the support of the angular measure. In practice, a sparse representation of the extremal dependence structure is obtained when only a few such groups of variables can be exhibited (compared to $2^{d}-1$) and/or when these groups involve a small number of variables (with respect to $d$). Here we develop this framework further, in order to propose a (soft) clustering technique in the region of extremes and derive effective $2$-d visual displays, sheding light on the structure of anomalies/extremes in sparse situations. This is achieved by modelling the distribution of extremes as a specific mixture model, where each component generates a different type $\alpha$ of extremes. In this respect, the present paper may be seen as an extension of Boldi and Davison, (2007); Sabourin and Naveau, (2014), where a Bayesian inference framework is designed for moderate dimensions ($d\leq 10$ say) and situations where the sole group of variables with the potential of being simultaneously large is $\{1,\ldots,d\}$ itself. In the context of mixture modelling (see e.g. Fruhwirth-Schnatter et al., (2018)), the Expectation-Maximization algorithm (EM in abbreviated form) permits to partition/cluster the set of extremal data through the statistical recovery of latent observations, as well as posterior probability distributions (inducing a soft clustering of the data in a straighforward manner) and, as a by-product, a similarity measure on the set of extremes: the higher the probability that their latent variables are equal, the more similar two extreme observations $X$ and $X^{\prime}$ are considered. The similarity matrix thus obtained naturally defines a weighted graph, whose vertices are the anomalies/extremes observed, paving the way for the use of powerful graph-mining techniques for community detection and visualization, see e.g. Schaeffer, (2007), Hu and Shi, (2015) and the references therein. Beyond its detailed description, the methodology proposed is applied to a real fleet monitoring dataset in the aeronautics domain and shown to provide useful tools for analyzing and interpreting abnormal data.

The paper is structured as follows. Basic concepts of MEV theory are briefly recalled in Section 2, in particular the concept of angular measure, together with the technique proposed in Goix et al., (2016, 2017) for estimating the (hopefully sparse) support of the latter, which determines the dependence structure of extremes arising from a heavy-tailed distribution. Section 3 details the mixture model we propose to describe the distribution of extreme data, based on the output of the support estimation procedure, together with the EM algorithm variant we introduce in order to estimate its parameters. It is next explained in Section 5 how to exploit the results of this inference method to define a similarity matrix of the extremal data, reflecting a weighted graph structure of the observed anomalies, and apply dedicated community detection and visualization techniques so as to extract meaningful information from the set of extreme observations. The relevance of the approach we promote is finally illustrated by numerical experiments, on synthetic and real data in Section 6. An implementation of the proposed method and the code for the experiments carried out in this paper are available online111 https://github.com/mchiapino/mevt_anomaly. Technical details are deferred to the Appendix section.

2 Background and Preliminaries

We start with recalling key notions of MEVT, the concept of angular measure in particular, as well as the inference method investigated in Goix et al., (2016, 2017) to estimate its support. Here and throughout, the Dirac mass at any point $x$ is denoted by $\delta_{x}$, the indicator function of any event $A$ by $\mathbbm{1}\{A\}$, the cardinality of any finite set $E$ by $|E|$. Capital letters generally refer to random quantities whereas lower case ones denote deterministic values. Finally, boldface letters denote vectors as opposed to Roman letters denoting real numbers.

2.1 Heavy-Tail Phenomena - Multivariate Regular Variation

Extreme Value Theory (EVT) describes phenomena that are not governed by an ’averaging effect’ but can be instead significantly impacted by very large values. By focusing on large quantiles rather than central statistics such as the median or the sample mean, EVT provides models for the unusual rather than the usual and permits to assess the probability of occurence of rare (extreme) events. Application domains are numerous and diverse, including any field related to risk management as finance, insurance, environmental sciences or aeronautics. Risk monitoring is a typical use case of EVT. The reader is referred to Coles, (2001) and the references therein for an introduction to EVT and its applications. In the univariate setting, typical quantities of interest are high quantiles of a random variable $X$, i.e. $1-p$ quantiles for $p\to 0$. When $p$ is of the same order of magnitude as $1/N$ or smaller, empirical estimates become meaningless. Another issue is the estimation of the probability of an excess over a high threshold $u$, $p_{u}=\mathbb{P}\left(\right)(X>u)$ when few (or none) observations are available above $u$. In such contexts, EVT essentially consists in using a parametric model (the generalized Pareto distributions) for the tail distribution, which is theoretically justified asymptotically, i.e. when $p\to 0$ or $u\to\infty$. Here and throughout we place ourselves in the context where the variable of interest is regularly varying; see Resnick, (1987, 2007) for a general introduction to regular variation and its applications to data analysis. In the univariate case the required assumption is the existence of a sequence $a_{n}>0$ such that $a_{n}\to\infty$ and a function $h(x)$ such that $n\mathbb{P}\left(X/a_{n}>x\right)\;\xrightarrow[n\rightarrow\infty]{}\;h(x)$, $x>0$. Notice that this assumption is satisfied by most textbook heavy tailed distributions, e.g. Cauchy, Student. In such a case $h$ is necessarily of the form $h(x)=Cx^{-\alpha}$ for some $C,\alpha>0$, where $\alpha$ is called the tail index of $X$ and $a_{n}$ may be chosen as $a_{n}=n^{1/\alpha}$. In the multivariate setting, consider a $d$-dimensional random vector $\mathbf{X}=(X_{1},\ldots,X_{d})$, the goal is to infer quantities such as $\mathbb{P}\left(X_{1}>x_{1},\ldots,X_{d}>x_{d}\right)$ for large $x_{1},\ldots,x_{d}$. A natural first step is to standardize each marginal distribution so that the $X_{j}$’s are all regularly varying with tail index $\alpha=1$ and scaling constant $C=1$. One convenient choice is to use the probability integral transform. For $\mathbf{x}=(x_{1},\ldots,x_{d})$, let $F_{j}(x_{j})=\mathbb{P}\left(X_{j}\leq x_{j}\right)$. Assuming that $F_{j}$ is continuous, the transformed variable $V_{j}=(1-F_{j}(X_{j}))^{-1}$ follows a Pareto distribution, $\mathbb{P}\left(V_{j}>v\right)=v^{-1}$, $v\geq 1$. In practice $F_{j}$ is unknown but its empirical version $\hat{F}_{j}$ may be used instead. Another option when each $X_{j}$ is regularly varying with tail index $\alpha_{j}$ is to estimate $(\alpha_{j},C_{j})$ using e.g. a generalized Pareto model above large thresholds, see Coles, (2001) or Beirlant et al., (2004) and the references therein. Then $V_{j}=X^{\alpha_{j}}/C_{j}$ is standard regularly varying, meaning that $n\mathbb{P}\left(V_{j}/n>x\right)\rightarrow x^{-1},x>0$ as $n\rightarrow\infty$. The multivariate extension of the latter assumption is that the standardized vector $\mathbf{V}=(V_{1},\ldots,V_{d})$ is regularly varying with tail index equal to $1$, i.e. there exists a limit Radon measure $\mu$ on $\mathbb{R}_{+}^{d}\setminus\{0\}$ such that

[TABLE]

for all $A$ in the continuity set of $\mu$ such that $0\notin\partial A$. The measure $\mu$ is called the exponent measure. In the standard setting characterized by (1), it is homogeneous of order $-1$, that is $\mu(tA)=t^{-1}\mu(A)$, where $tA=\{t\mathbf{v},\mathbf{v}\in A\},A\subset\mathbb{R}_{d}^{+}$ and $\mu\{\mathbf{x}\in\mathbb{R}_{+}^{d}:x_{j}\geq 1\}=1$. Assumption (1) applies immediately to the problem of estimating the probability of reaching a set $tA$ which is far form ${\mathbf{0}}$ (i.e. $t$ is large): one may write $\mathbb{P}\left(\right)(\mathbf{V}\in tA)\approx\frac{1}{t}\mu(A)$, so that estimates of $\mu$ automatically provide estimates for such quantities. In a word, $\mu$ may be used to characterize the distributional tail of $\mathbf{V}$. For modeling purposes, the homogeneity property $\mu(t\cdot)=t^{-1}\mu(\cdot)$ suggests a preliminary decomposition of $\mu$ within a (pseudo)-polar coordinates system, as detailed next.

2.2 Angular Measure - Dependence in the Extremes

Consider the sum-norm $\|{\mathbf{v}}\|:=v_{1}+\ldots+v_{d}$ and $\mathcal{S}_{d}:=\{\mathbf{w}\in\mathbb{R}_{+}^{d}:\|{\mathbf{w}}\|=1\}$ the $d$-dimensional simplex. Introduce the polar transformation $T:\mathbf{v}\mapsto T(\mathbf{v})=(r,\mathbf{w})$ defined on $\mathbb{R}_{+}^{d}\setminus\{\mathbf{0}\}$, where $r=\|\mathbf{v}\|$ is the radial component and $\mathbf{w}={r}^{-1}{\mathbf{v}}$ is the angular one. Now define the angular measure $\Phi$ on $\mathcal{S}_{d}$ (see e.g. Resnick, (2007) or Beirlant et al., (2004) and the references therein): $\Phi(A):=\mu\left\{\mathbf{v}:\|{\mathbf{v}}\|>1,\|{\mathbf{v}}\|^{-1}\mathbf{v}\in A\}\right\}$, with $A\subset\mathcal{S}_{d}$. Notice that $\Phi(\mathcal{S}_{d})<\infty$ and, by homogeneity,

[TABLE]

In other words the exponent measure $\mu$ factorizes into a tensor product of a radial component and an angular component. Setting $R=\|\mathbf{V}\|$ and $\mathbf{W}=R^{-1}\mathbf{V}$, a consequence is that

[TABLE]

for all measurable set $A\subset\mathcal{S}_{d}$ such that $\Phi(\partial A)=0$ and $r>1$. Hence, given that the radius $R$ is large, $R$ and the angle $\mathbf{W}$ are approximately independent, the distribution of $\mathbf{W}$ is approximately the angular measure – up to a normalizing constant $\Phi(\mathcal{S}_{d})$ – and $R$ follows approximately a Pareto distribution. As it describes the distribution of the directions formed by the largest observations, the angular measure exhaustively accounts for the dependence structure in the extremes. Our choice of a standard regular variation framework (1) and that of the sum-norm yield the following moment constraint on $\Phi$:

[TABLE]

In addition, the normalizing constant is explicit:

[TABLE]

Remark 1**.**

The choice of the sum-norm here is somewhat arbitrary. Any other norm on $\mathbb{R}^{d}$ for the pseudo-polar transformation is equally possible, leading to alternative moment constraints and normalizing constants. The advantage of the sum-norm is that it allows convenient probabilistic modeling of the angular component $\mathbf{w}$ on the unit simplex.

2.3 Support Estimation - the DAMEX Algorithm

We now expose the connection between $\Phi$’s (or equivalently, $\mu$’s) support and the subsets of components which may simultaneously take very large values, while the others remain small.

Sparse support. Fix $\alpha\subset\{1,\ldots,d\}$ and consider the associated truncated cone

[TABLE]

The family $\{\mathcal{C}_{\alpha},\alpha\subset\{1,\ldots,d\},\alpha\neq\emptyset\}$ first introduced in Goix et al., (2016) defines a partition of $\mathbb{R}_{+}^{d}\setminus[0,1]^{d}$ which is of particular interest for our purpose: notice first that, by homogeneity of $\mu$, the following equivalence holds: $\Phi(\mathcal{S}_{\alpha})>0\Leftrightarrow\mu(\mathcal{C}_{\alpha})>0$, where $\mathcal{S}_{\alpha}=\left\{\mathbf{v}\in\mathbb{R}_{+}^{d}:\;\|{v}\|=1,\;v_{i}>0\text{ for }i\in\alpha,v_{i}=0\text{ for }i\notin\alpha\right\}$, $\varnothing\neq\alpha\subset\{1,\ldots,d\}$. Observe next that $\mu(\mathcal{C}_{\alpha})>0$ means that the limiting rescaled probability that ‘all features in $\alpha$ are simultaneously large, while the others are small’ is non zero. Precisely, consider the $\epsilon$-thickened rectangle

[TABLE]

which corresponds to the event that all features in $\alpha$ are large, while the other are small. The $\mathcal{R}_{\alpha}^{\epsilon}$’s define again a partition of $\mathbb{R}_{+}^{d}\setminus[0,1]^{d}$ for each fixed $\epsilon\geq 0$. In addition, we have that $\mathcal{C}_{\alpha}=\cap_{\epsilon>0,\epsilon\in\mathbb{Q}}\mathcal{R}_{\alpha}^{\epsilon}$, so that by upper continuity of $\mu$,

[TABLE]

with

[TABLE]

In the sequel, set $\mu_{\alpha}=\mu(\mathcal{C}_{\alpha})$, \mathbb{M}=\big{\{}\alpha\subset\{1,\ldots,d\},\alpha\neq\varnothing,\mu_{\alpha}>0\big{\}}. Although every $\mu_{\alpha}$ may be positive in theory, a reasonable assumption in many practical high dimensional situations is that $\mu_{\alpha}=0$ for the vast majority of the $2^{d}-1$ cones $\mathcal{C}_{\alpha}$. In other words, not all combinations of coordinates of $\mathbf{V}$ can be large together, so that the support of $\mu$ (and that of $\Phi$) is sparse.

Support estimation. The task of estimating the support of $\mu$ (or $\Phi$) has recently received increasing attention in the statistics and machine learning literature. Chautru, (2015) first proposed a non parametric clustering approach involving principal nested spheres, which provides great flexibility at the price of computational cost. In contrast, Goix et al., (2016)’s methods rely on the above mentioned partition of the unit sphere into $2^{d}-1$ sub-simplices $\Omega_{\alpha}$ and led to the so-called DAMEX algorithm which computational complexity $O(dn\log n)$ scales well with higher dimensions. Their algorithm produces the list of $\alpha$’s such that the empirical counterpart of $\mu_{\alpha}$ (denoted $\hat{\mu}_{\alpha}$ in the sequel) is non zero. Defining a threshold $m_{\min}>0$ below which $\hat{\mu}_{\alpha}$ is deemed as negligible, one thus obtains a list of subsets $\widehat{\mathbb{M}}=\{\alpha\subset\{1,\ldots,d\}:\hat{\mu}_{\alpha}>\mu_{\min}\}$. A uniform boud on the error $|\hat{\mu}_{\alpha}-\mu_{\alpha}|$ is derived in Goix et al., (2017) which scales roughly as $k^{-1/2}$, where $k$ is the order of magnitude of the number of largest observations used to learn $\mathbb{M}$ and the $\mu_{\alpha}$’s. In Simpson et al., (2018) the original DAMEX framework is refined in order to also model extremes in the directions $\Omega_{\alpha}$ where the angular measure does not concentrate. A third algorithm named CLEF has been proposed by Chiapino and Sabourin, (2016) which allows to cluster together different sub-simplices $\Omega_{\alpha}$’s which are close in terms of symmetric difference of the subsets $\alpha$’s. This is particularly useful in situations where the empirical angular mass is scattered onto a large number of sub-simplices, so that DAMEX fails to exhibit a list $\hat{\mathbb{M}}$ of reasonable size. Asymptotic guarantees for the latter approach and variants, leading to statistical tests with controllable asymptotic type I error are derived in Chiapino et al., (2018).

In the present paper, support estimation is only a preliminary step before mixture modeling. We decided to use DAMEX in view of its computational simplicity and the statistical guarantees it offers, considering the fact that its output was very similar to CLEF’s on the aeronautics dataset considered in our usecase (see Section 6.2). Using the above mentioned alternatives is certainly possible but for the sake of brevity we shall only present the results obtained using DAMEX as a preprocessing step. We now briefly describe how DAMEX works.

**The DAMEX algorithm. ** Given a dataset $(\mathbf{X}_{i})_{i\leq n}$ of independent data distributed as $\mathbf{X}$, DAMEX proceeds as follows. First, replace the unknown marginal distributions $F_{j}$ with their empirical counterpart $\hat{F}_{j}(x)=\frac{1}{n}\sum\mathbbm{1}\{X_{i,j}<x\}$ and define next $\hat{V}_{i,j}=(1-\hat{F}_{j}(X_{i,j}))^{-1}$ and $\mathbf{\hat{V}}_{i}=(\hat{V}_{i,1},\ldots,\hat{V}_{i,d})$. Then choose some $k\ll n$ large enough (typically $k=O(\sqrt{n})$) and define $\hat{\mu}_{\alpha}$ as the empirical counterpart of $\mu(R_{\alpha}^{\epsilon})$ with $t$ replaced by $n/k$, that is $\hat{\mu}_{\alpha}=(1/k)\sum_{i=1}^{n}\mathbbm{1}\{\hat{\mathbf{V}_{i}}\in\frac{n}{k}\mathcal{R}_{\alpha}^{\epsilon}\}$. Notice that the above description is a variant of the original algorithm in Goix et al., (2016) which uses thickened cones $\mathcal{C}_{\alpha}^{\epsilon}$ instead of $\mathcal{R}_{\alpha}^{\epsilon}$. However finite sample guarantees in Goix et al., (2017) are obtained using the latter rather than the original $\mathcal{C}_{\alpha}^{\epsilon}$’s, that is why using the $\mathcal{R}_{\alpha}^{\epsilon}$’s is preferred.

3 A Mixture Model For Multivariate Extreme Values

The purpose of this section is to develop a novel mixture model for the angular distribution $\Phi$ of the largest instances of the dataset, indexed by $\alpha\in\mathbb{M}$, where $\mathbb{M}$ is $\Phi$’s support. Each component $\alpha\in\mathbb{M}$ of the mixture generates instances $\mathbf{V}$ such that $V_{j}$ is likely to be large for $j\in\alpha$ and the latent variables of the model take their values in $\mathbb{M}$. In practice, we adopt a plug-in approach and identify $\mathbb{M}$ with $\widehat{\mathbb{M}}$, the output of DAMEX. As the distribution of extremes may be entirely characterized by the distribution of their angular component $\mathbf{W}\in\mathcal{S}_{d}$ (see the polar decompositions (2) and (3)), a natural model choice is that of Dirichlet mixtures. We next show how to design a ’noisy’ version of the model for subasymptotic observations and how to infer it by means of an EM procedure based on a truncated version of the original dataset, surmounting difficulties related to the geometry of $\Phi$’s support.

3.1 Angular Mixture Model for the Directions along which Anomalies Occur

Recall from Section 2.3 that $\mathcal{S}_{d}$ is naturally partitioned into $2^{d}-1$ sub-simplices $\mathcal{S}_{\alpha}$. Our key assumption is that the support of $\mu$ (or $\Phi$) is sparse in the sense that $|\mathbb{M}|\ll 2^{d}$, where $\mathbb{M}=\{\alpha:\mu(\mathcal{C}_{\alpha})>0\}=\{\alpha:\Phi(\mathcal{S}_{\alpha})>0\}$. Let $K$ denote the number of subsets $\alpha\in\mathbb{M}$ of cardinality at least $2$ and let $d_{1}\in\{0,\ldots,d\}$ be the number of singletons $\{j\}\in\mathbb{M}$. Without loss of generality we assume that these singletons correspond to the first $d_{1}$ coordinates, so that $\mathbb{M}=\{\alpha_{1},\ldots,\alpha_{K},\{1\},\ldots,\{d_{1}\}\}$. For simplicity, we also suppose that the sets $\alpha\in\mathbb{M}$ are not nested, an hypothesis which can be relaxed at the price of additional notational complexity. In view of (5), the angular measure then admits the decomposition

[TABLE]

where $\Phi_{\alpha_{k}}$ is a probability measure on $\mathcal{S}_{\alpha_{k}}$, the weights $\pi_{k}$ satisfy $\sum_{k\leq K+j}\pi_{k}=1$ and $\mathbf{e}_{j}=(0,\ldots,1,\ldots,0)$ is the $j^{th}$ canonical basis vector of $\mathbb{R}^{d}$. The singletons weights derive immediately from the moment constraint (4): for $i\leq d_{1}$,

[TABLE]

We obtain

[TABLE]

where the vector $\mathbf{\pi}\in[0,1]^{K+d_{1}}$ must satisfy

[TABLE]

Equation (7) determines the structure of the angular distribution of the largest observations. For likelihood-based inference, a parametric model for each component $\Phi_{\alpha_{k}}$ of the angular measure must be specified. One natural model for probability distributions on a simplex is the Dirichlet family, which provides a widely used prior in Bayesian statistics for data clustering purposes in particular. We recall that the Dirichlet distribution on a simplex $\mathcal{S}_{\alpha}$ admits a density $\varphi_{\alpha}$ with respect to the $(|\alpha|-1)$-dimensional Lebesgue measure which is denoted by $\mathrm{d}\mathbf{w}$ for simplicity. It can be parameterized by a mean vector $\mathbf{m}_{\alpha}\in\mathcal{S}_{\alpha}$ and a concentration parameter $\nu_{\alpha}>0$, so that for $\mathbf{w}\in\mathcal{S}_{\alpha}$,

[TABLE]

Refer to e.g. Müller and Quintana, (2004) for an account of Dirichlet processes and mixtures of Dirichlet Processes applied to Bayesian nonparametrics. We emphasize that our context is quite different: a Dirichlet Mixture is used here as a model for the angular component of the largest observations, not as a prior on parameters. This modeling strategy for extreme values was first proposed in Boldi and Davison, (2007) and revisited in Sabourin and Naveau, (2014) to handle the moment constraint (4) via a model re-parametrization. In both cases, the focus was on moderate dimensions. In particular, both cited references worked under the assumption that the angular measure concentrates on the central simplex $\Omega_{\{1,\ldots,d\}}$ only. In this low dimensional context, the main purpose of the cited authors was to derive the posterior predictive angular distribution in a Bayesian framework, using a variable number of mixture components concentrating on $\Omega_{\{1,\ldots,d\}}$. Since the set of Dirichlet mixture distributions with an arbitrary number of components is dense among all probability densities on the simplex, this model permits in theory to approach any angular measure for extremes. The scope of the present paper is different. Indeed we are concerned with high dimensional data (say $d\simeq 100)$ and consequently we do not attempt to model the finest details of the angular measure. Instead we intend to design a model accounting only for information which is relevant for clustering. Since an intuitive summary of an extreme event in a high dimensional context is the subset $\alpha$ of features it involves, we assign one mixture component per sub-simplex $\Omega_{\alpha}$ such that $\alpha\in\mathbb{M}$. Thus we model each $\Phi_{\alpha}$ by a single Dirichlet distribution with unknown parameters $m_{\alpha},\nu_{\alpha}$. Using the standard fact that for such a distribution, $\int_{\mathbb{S}_{\alpha}}\mathbf{w}\varphi_{\alpha}(\mathbf{w}|m_{\alpha},\nu_{\alpha})\mathrm{d}\mathbf{w}=\mathbf{m}_{\alpha}$, the moment constraint (4) becomes:

[TABLE]

where $\mathbf{m}_{k}=\mathbf{m}_{\alpha_{k}}$ for $k\leq K$.

3.2 A Statistical Model for Large but Sub-asymptotic Observations.

Recall from (3) that $\Phi$ is the limiting distribution of $\mathbf{V}$ for large $R$’s. In practice, we dispose of no realization of this limit probability measure and the observed angles corresponding to radii $R>r_{0}$ follow a sub-asymptotic version of $\Phi$. In particular, if the margins $V_{j}$ have a continuous distribution, we have $\mathbb{P}\left(\right)(V_{j}\neq 0)=1$ so that all the $\mathbf{V}_{i}=(V_{i,1},\ldots,V_{i,d})$, $1\leq i\leq n$, lie in the central cone $\mathcal{C}_{\{1,\ldots,d\}}$ (this is also true using the empirical versions $\hat{\mathbf{V}}_{i}$ defined in subsection 2.3). In the approach we propose, the deviation of $\mathbf{V}$ from its asymptotic support, which is $\bigcup_{\alpha\in\mathbb{M}}\mathcal{C}_{\alpha}$, is accounted for by a noise $\bm{\varepsilon}$ with light tailed distribution, namely an exponential distribution. That is, we assume that $\mathbf{V}=R\,\mathbf{W}+\bm{\varepsilon}$, see Model 1 below. As is usual for mixture modeling purposes, we introduce a multinomial latent variable $\mathbf{Z}=(Z_{1},\;\ldots,\;Z_{K+d_{1}})$ such that $\sum_{k}Z_{k}=1$ and $Z_{k}=1$ if $\mathbf{W}$ has been generated by the $k^{th}$ component of the angular mixture (7). In a nutshell, the type of anomaly/extreme is encoded by the latent vector $\mathbf{Z}$. Then, for $k\leq K$, $\mathbb{P}\left(Z_{k}=1\right)=\pi_{k}$, while, for $K<k\leq K+d_{1}$, $\mathbb{P}\left(Z_{k}=1\right)=d^{-1}$. The unknown parameters of the model are $\bm{\theta}=(\mathbf{\pi},\mathbf{m},\bm{\nu})$, where $\nu_{k}>0$ and $\mathbf{\pi}=(\pi_{1},\ldots,\pi_{K}),\;\mathbf{m}=(\mathbf{m}_{1},\ldots,\mathbf{m}_{K})$ must satisfy the constraints (8) and (9), as well as the exponential rates $\bm{\lambda}=(\lambda_{1},\;\ldots,\;\lambda_{K+d_{1}})$, where $\lambda_{k}>0$. Figure 1 illustrates Model 1 in dimension $d=3$.

Model 1** (Sub-asymptotic mixture model).**

•

Consider a standard regularly varying random vector ${\mathbf{V}}$ satisfying (1) (typically $V_{j}=(1-\hat{F}_{j}(X_{j}))$ for $\hat{F}_{j}$ an estimate of the marginal distribution $F_{j}$ of $X_{j}$, see subsection 2.1).

•

Let $R=\|{{\mathbf{V}}}\|$. Fix some high radial threshold $r_{0}$, typically a large quantile of the observed radii. Let $\mathbf{Z}$ be a hidden variable indicating the mixture component in (7). Conditionally to $\{R>r_{0},Z_{k}=1\}$, ${\mathbf{V}}$ decomposes as

$\displaystyle{\mathbf{V}}\;=\;\mathbf{V}_{k}+\bm{\varepsilon}_{k}\;=\;R_{k}\mathbf{W}_{k}+\bm{\varepsilon}_{k},$

(10)

where $\mathbf{V}_{k}\in\mathcal{C}_{\alpha_{k}}$, $\bm{\varepsilon}_{k}\in\mathcal{C}_{\alpha_{k}}^{\perp}$, $R_{k}=\|\mathbf{V}_{k}\|$, $\mathbf{W}_{k}=R_{k}^{-1}\mathbf{V}_{k}\in\mathcal{S}_{\alpha_{k}}$. The components $R_{k},\mathbf{W}_{k},\bm{\varepsilon}_{k}$ are independent from each other. The noise’s components are i.i.d. according to a translated exponential distribution with rate $\lambda_{k}$, $R_{k}$ is Pareto distributed above $r_{0}$ and $\mathbf{W}_{k}$ is distributed as $\Phi_{k}$, that is

$\left\{\begin{aligned} &\mathbb{P}\left(R_{k}>r\right)=r_{0}r^{-1},r>r_{0}\,,\\ &\mathbf{W}_{k}\sim\Phi_{k}\,,\\ &\varepsilon_{j}\sim 1+\mathcal{E}\mathit{xp}(\lambda_{k}),j\in\{1,\ldots,d\}\setminus\alpha_{k}\,,\end{aligned}\right.$

with $\Phi_{k}=\varphi_{k}(\,\cdot\,|m_{k},\nu_{k})$ if $k\leq K$, and $\Phi_{k}=\delta_{\mathbf{e}_{k-K}}$ if $K<k\leq K+d_{1}$.

4 Statistical Inference via EM Algorithm.

In the mixture model setting described above with hidden variables $Z_{i}$, likelihood optimization is classically performed using an EM algorithm (Dempster et al.,, 1977). This method consists in performing in turn the so-called E-step and M-step at each iteration $t$. Denoting by $\bm{\theta}_{t}$ the value at iteration $t$ of the set of unknown model parameters, the posterior probabilities

[TABLE]

are computed during the E-step and define the objective function

[TABLE]

The latter serves as a proxy for the log-likelihood and can be maximized with respect to $\bm{\theta}$ with standard optimization routines during the M-step, which yields $\bm{\theta}_{t+1}=\operatorname*{arg\,max}_{\bm{\theta}}Q(\bm{\theta},\gamma^{(t)})$. The procedure stops when the value $Q(\bm{\theta}_{t},\gamma^{(t)})$ reaches a stationary point and the latest pair $(\bm{\theta}_{t},\gamma^{(t)})$ is returned.

The likelihood for Model 1, $p({\mathbf{v}}|\bm{\theta}=(\mathbf{m},{\bm{\nu}},\mathbf{\pi},\bm{\lambda}))$, for one observation ${\mathbf{v}}\in(1,\infty)^{d},\|{{\mathbf{v}}}\|\geq r_{0}$, follows directly from the model specification,

[TABLE]

where $f_{\varepsilon}(\,\cdot\,\;|\lambda_{k})$ denotes the marginal density of the noise $\bm{\varepsilon}_{k}$ given the noise parameter $\lambda_{k}$. As specified in Model 1, in this paper we set $f_{\varepsilon}(x|\lambda_{k})=\lambda_{k}e^{-\lambda_{k}(x-1)}$, $x>1$ (a translated exponential density), but any other light tailed distribution could be used instead. Notice that the term $r_{k}^{-|\alpha_{k}|-1}=r_{k}^{-2}r_{k}^{-|\alpha_{k}|+1}$ is the product of the radial Pareto density and the Jacobian term for the change of variables $T_{k}:\mathbf{V}_{k}\mapsto(R_{k},\mathbf{W}_{k})$. Recall that the constraints are

[TABLE]

and that $\mathbf{\pi}=(\pi_{1},\ldots,\pi_{K})$ and $\mathbf{m}=(\mathbf{m}_{1},\ldots,\mathbf{m}_{K})$ satisfy (8) and (9). The latter linear constraint on $(\bm{\pi},\mathbf{m})$ implies that $\mathbf{m}$ and $\bm{\pi}$ cannot be optimized independently, which complicates the M-step of an EM-algorithm. Thus we begin with a re-parametrization of the model ensuring that the moment constraint (4) is automatically satisfied.

Re-parametrization. In a lower dimensional Bayesian framework, earlier works (Sabourin and Naveau, (2014)) have proposed a re-parametrization of the pair $(\bm{\pi},\bf m)$ ensuring that the moment constraint (4) is automatically satisfied. This consists in a sequential definition of the mixture centers $m_{k}$ together with the involving partial barycenters of the remaining components $(m_{k+1},\ldots,m_{k})$. The advantage if this construction is that the resulting parameter has a intuitive interpretation which facilitates the definition of a prior, while allowing for efficient MCMC with reversible jumps sampling (Green, (1995)) of the posterior distribution. However, how to adapt this re-parameterization to our context where several sub-simplices are involved remains an open question and we did not pursue this idea further. The re-parametrization that we propose here consists in working with the product parameter $\rho_{k,j}=\pi_{k}m_{k,j}$ instead of the pair $(\pi_{k},m_{k,j})$. Namely, consider a $K\times(d-d_{1})$ matrix $\bm{\rho}=({\bm{\rho}_{1}}^{\top},\ldots,{\bm{\rho}_{K}}^{\top})$ where $\rho_{k,j}>0$ for $j\in\alpha_{k}$ and $\rho_{k,j}=0$ otherwise. Then, for all $k\in\{1\ldots,K\}$, set

[TABLE]

Then (8) and (9) together are equivalent to

[TABLE]

In the sequel we denote respectively by $p({\mathbf{v}}|\bm{\rho},\bm{\nu},\bm{\lambda}):=p({\mathbf{v}}|\bm{\pi},\mathbf{m},\bm{\nu},\bm{\lambda})$ and $\varphi_{k}(\mathbf{w}|\bm{\rho}_{k},\nu_{k}):=\varphi_{k}(\mathbf{w}|\mathbf{m}_{k},\nu_{k})$ the likelihood and the Dirichlet densities in the re-parameterized model, where $(\mathbf{m},\bm{\pi})$ are obtained from $\bm{\rho}$ via (13). By abuse of notations, let $\bm{\theta}$ denote in the sequel the set of parameters of the re-parameterized version of Model 1, that is $\bm{\theta}=(\bm{\rho},\bm{\nu},\bm{\lambda})$, and let $\Theta$ be the parameter space, that is the set of $\bm{\theta}$’s such that constraints (12) and (14) hold.

EM algorithm. We summarize below the EM algorithm in our framework. Let $n_{0}\leq n$ be the number of observations ${\mathbf{V}}_{i}$ such that $\|{{\mathbf{V}}_{i}}\|>r_{0}$. To alleviate notations, we may relabel the indices $i$ so that these observations are ${\mathbf{V}}_{1:n_{0}}=({\mathbf{V}}_{1},\ldots,{\mathbf{V}}_{n_{0}})$. Let $\mathbf{Z}_{i}=(Z_{i,1},\ldots,Z_{i,K+d_{1}}),i\leq n_{0}$ be the hidden variables associated with ${\mathbf{V}}_{1:n_{0}}$. Also let $p({\mathbf{v}}|\bm{\theta},z_{k}=1)$ denote the conditional density of ${\mathbf{V}}$ given $(Z_{k}=1,\bm{\theta})$. In view of the likelihood (11), it is given by

[TABLE]

**EM algorithm for Model 1 **

Input Extreme standardized data ${\mathbf{V}}_{1:n_{0}}$.

•

Initialization Choose a starting value for $\bm{\theta}$ (See Remark 2).

•

Repeat until convergence:

E-step: compute for $1\leq i\leq n_{0}$ and $k\leq K+d_{1}$, $\gamma_{i,k}=\mathbb{P}\left(Z_{i,k}=1\;\middle|\;{}{\mathbf{V}}_{i},\bm{\theta}\right)$ according to (17). Set $\bm{\gamma}=(\gamma_{i,k})_{i\leq n_{0},k\leq K+d_{1}}$.

M-step: Solve the optimization problem $\max_{\bm{\theta}\in\Theta}Q(\bm{\theta},\bm{\gamma})$where Q(\bm{\theta},\bm{\gamma})=\sum_{i=1}^{n_{0}}\sum_{k=1}^{K+d_{1}}\gamma_{i,k}\big{(}\log\pi_{k}+\log p({\mathbf{V}}_{i}|\bm{\theta},z_{i,k}=1)\big{)} is a lower bound for the likelihood and $\pi_{k}=\mathbb{P}\left(\right)(Z_{i,k}=1|\bm{\theta})$, i.e.

$\pi_{k}=\begin{cases}\sum_{\ell\in\alpha_{k}}\rho_{k,l}&\text{ for }1\leq k\leq K\;,\\ d^{-1}&\text{ for }K<k\leq K+d_{1}\;,\end{cases}$

(16)

where $p({\mathbf{V}}_{i}|\bm{\theta},z_{i,k}=1)$ is given by (15).Denote by $\bm{\theta}^{\star}$ the solution, set $\bm{\theta}=\bm{\theta}^{\star}$.

Remark 2**.**

*In this work the starting values for the concentration parametrers $\nu_{k}$ are set to $20$, those for the exponential rates are set to $\lambda_{k}=0.01$. Finally, one may easily construct a matrix $\bm{\rho}$ satisfying the constraint (14) starting with any matrix $\tilde{\bm{\rho}}$ such that $\tilde{\rho}_{k,j}=0$ for $j\notin\alpha_{k}$ and $\tilde{\rho}_{k,j}>0$ otherwise, and then defining $\bm{\rho}$ via $\rho_{k,j}=(\sum_{l=1}^{K}\tilde{\rho}_{l,j})^{-1}\tilde{\rho}_{k,j}$. *

We now describe at length the E-step and the M-step of the algorithm.

E-step. The $\gamma_{i,k}$’s are obtained using the Bayes formula, for $1\leq k\leq K+d_{1}$,

[TABLE]

where $\pi_{k}$ is defined in (16) and $p({\mathbf{V}_{i}}|Z_{i,k}=1,\bm{\theta})$ is given by (15).

M-step. Here optimization of $Q(\bm{\theta},\gamma)$ with respect ot $\bm{\theta}=(\bm{\rho},\bm{\nu},\bm{\lambda})$ is performed under constraints (12), (14). Since $Q$ decomposes into a function of $(\bm{\rho},\bm{\nu})$ and a function of $\bm{\lambda}$, and since the constraints on $\bm{\rho},\bm{\nu}$ and $\bm{\lambda}$ are independent, maximization can be performed independently over the two blocks. Indeed, gathering terms not depending on $\bm{\theta}$ into a constant $C$,

[TABLE]

where

[TABLE]

Here we set $\alpha_{k}=\{k-K\}$ for $K<k\leq K+d_{1}$, in accordance with the notations from Section 3.1. Notice that the dependence of $Q_{1}$ and $Q_{2}$ on $\bm{\gamma}$ is omitted for the sake of concision. With these notations

[TABLE]

The function $Q_{1}$ being non-concave we use the python package mystic (McKerns et al., (2012)) to maximize it. For our choice of translated exponential noise, $f_{\varepsilon}(v|\lambda_{k})=\lambda_{k}e^{-\lambda_{k}(v-1)}$, $v\geq 1$, the maximizer of $Q_{2}$ has an explicit expression,

[TABLE]

Remark 3**.**

Let $\bm{\gamma}^{t}$ and $\bm{\theta}^{t}$ be the results of the $t$-th iteration of the algorithm then we conclude the iterative process if $Q(\bm{\theta}^{t},\bm{\gamma}^{t})<Q(\bm{\theta}^{t-1},\bm{\gamma}^{t-1})+\epsilon$, with $\epsilon$ a small threshold.

5 Graph-based Clustering and Visualization Tools

Beyond the hard clustering that may be straightforwardly deduced from the computation of the likeliest values $z_{1},\;\ldots,\;z_{n_{0}}$ for the hidden variables given the ${\mathbf{V}}_{i}$’s and the parameter estimates produced by the EM algorithm, the statistical model previously introduced defines a natural structure of undirected weighted graph on the set of observed extremes, which interpretable layouts (graph drawing) can be directly derived using classical solutions. Indeed, a partition (hard clustering) of the set of (standardized) anomalies/extremes ${\mathbf{V}}_{1},\;\ldots,\;{\mathbf{V}}_{n_{0}}$ is obtained by assigning membership of each ${\mathbf{V}}_{i}$ in a cluster (or cone/sub-simplex ) determined by the component of the estimated mixture model from which it arises with highest probability: precisely, one then considers that the abnormal observation ${\mathbf{V}}_{i}$ is in the cluster indexed by

[TABLE]

and is of type $\alpha_{k_{i}}$. However, our model-based approach brings much more information and the vector of posterior probabilities $(\gamma_{i,1},\;\ldots,\;\gamma_{i,K+d_{1}})$ output by the algorithm actually defines soft membership and represent the uncertainty in whether anomaly ${\mathbf{V}}_{i}$ is in a certain cluster. It additionally induces a similarity measure between the anomalies: the higher the probability that two extreme values arise from the same component of the mixture model, the more similar they are considered. Hence, consider the undirected graph whose vertices, indexed by $i=1,\;\ldots,\;n_{0}$, correspond to the extremal observations ${\mathbf{V}}_{1},\;\ldots,\;{\mathbf{V}}_{n_{0}}$ and whose edgeweights are $w_{{\bm{\theta}}}({\mathbf{V}}_{i},{\mathbf{V}}_{j})$, $1\leq i\neq j\leq n_{0}$, where

[TABLE]

Based on this original graph description of the set of extremes, it is now possible to rank all anomalies (i.e. extreme points) by degree of similarity to a given anomaly ${\mathbf{V}}_{i}$

[TABLE]

and extract neighborhoods $\{{\mathbf{V}}_{(i,1)},\;\ldots,\;{\mathbf{V}}_{(i,l)}\}$, $l\leq n_{0}$.

Graph-theoretic clustering. We point out that many alternative methods to that consisting in assigning to each any anomaly/extreme its likeliest component (i.e. model-based clustering) can be implemented in order to partition the similarity graph thus defined into subgraphs whose vertices correspond to similar anomalies, ranging from tree-based clustering procedures to techniques based on local connectivity properties through spectral clustering. One may refer to e.g. Schaeffer, (2007) for an account of graph-theoretic clustering methods.

Graph visualization. In possible combination with clustering, graph visualization techniques (see e.g. Hu and Shi, (2015)), when the number $n_{0}$ of anomalies to be analyzed is large, can also be used to produce informative layouts. Discussing the merits and limitations of the wide variety of approaches documented in the literature in this purpose is beyond the scope of this paper. The usefulness of the weighted graph representation proposed above combined with state-of-the-art graph-mining tools is simply illustrated in Section 6.2 and 6.3. We point out however that alternatives to the (force-based) graph drawing method used therein can be naturally considered, re-using for instance the eigenvectors of the graph Laplacian computed through a preliminary spectral clustering procedure (see e.g. Athreya et al., (2017) and the references therein for more details on spectral layout methods).

6 Illustrative Experiments

The aim of our experiments is double. First, investigate the goodness of fit of the Dirichlet mixture model fitted via the EM algorithm on simulated data from the model. Second, provide empirical evidence of the relevance of the approach we promote for anomaly clustering/visualization with real world data. Comparisons with state-of-the-art methods standing as natural competitors are presented for this purpose.

6.1 Experiments on Simulated Data

To assess the performance of the proposed estimator of the dependence structure and of the EM algorithm, we generate synthetic data according to Model 1. The dimension is fixed to $d=100$ and the mixture components, that is the elements of $\mathbb{M}=\{\alpha_{1},\ldots,\alpha_{K}\}$, are randomly chosen in the power set of $\{1,\ldots,d\}$ with $K=50$. The coefficients of the matrix $\rho$ which determines the weights and centers through Eq. (14) in the Supplementary Material is also randomly chosen, then its columns are normalized so that the moment constraint is satisfied. Finally. we fix $\nu_{k}=20$ for $1\leq k\leq K$ and $\lambda_{k}$, $1\leq k\leq K+d_{1}$, are successively set to $1$, $0.75$, $0.5$, $0.25$ and $0.1$ to vary the noise level in the experiments. Then each point ${\mathbf{V}}_{i}=R_{i}\mathbf{W}_{i}+\bm{\varepsilon}_{i}$, $i\leq n$, is generated with probability $\pi_{k},k\in\{1,\ldots,K\}$ according to the mixture component $k\leq K$, that is

[TABLE]

and with probability $\frac{1}{d}$ according to component $k\in\{K,\ldots,K+d_{1}\}$ in such a way that

[TABLE]

The threshold $r_{0}$ above which points are considered as extreme is fixed to $100$. On this toy example, the pre-processing step that consists in applying DAMEX for recovering $\mathbb{M}$ produces an exact estimate, so that $\hat{\mathbb{M}}=\mathbb{M}$. Then the procedure described in Algorithm 4 is applied. Tables 1 and 2 show the average absolute errors for the estimates $\widehat{\rho}$, $\widehat{\nu}$ and $\widehat{\lambda}$ on $50$ datasets of the $n_{0}$ generated extreme points, for $n_{0}=1e+3,2e+3$, namely

[TABLE]

On this toy example, estimates of the means and weights, as well as those of the noise parameters, are almost exact. In contrast, the estimates of the $\nu_{k}$’s are not that accurate, but, as shown next, this drawback does not jeopardize cluster identification.

The performance in terms of cluster identification is measured as follows: for each point ${\bm{v}}_{i}$, the true label $y_{i}\in\{1,\ldots,K+d_{1}\}$ is compared with the label obtained via assignment to the highest probable component, that is $\hat{y}_{i}=\operatorname*{arg\,max}_{k\in\{1,\;\ldots,\;K+d_{1}\}}\gamma_{i,k}$. Table 3 shows the average number of labeling errors for different values of $n_{0}$ and $\lambda_{k}$.

Figure 2 illustrates the relevance of the proposed approach regarding anomaly visualization. A test set of size $100$ consisting of extreme data is simulated as above, and the corresponding matrix $\widehat{\bm{\gamma}}$ is computed according to (17) with ${\bm{\theta}}$ taken as the output of the training step (i.e. Algorithm 4 run with the training dataset of $n_{0}=2e3$ points). Finally an adjacency matrix $w_{\widehat{\bm{\theta}}}({\mathbf{v}}_{i},{\mathbf{v}}_{j})$ is obtained as detailed in Section 5, on which we apply spectral clustering in order to group the points according to the similarities measured by $w$. Graph visualization of $w$ is next performed using the python package ’Networkx’ Hagber et al., (2008), that provides a spring layout of the graph according to the Früchtermen-Reingold algorithm, see Früchterman and Reingold, (1991). A hard thresholding is applied to the edges in $w$ in order to improve readability: edges $(i,j)$ such that $w_{\widehat{\bm{\theta}}}({\mathbf{v}}_{i},{\mathbf{v}}_{j})<\epsilon$ with $\epsilon$ a small threshold are removed. Each cluster output by the spectral clustering method is identified with a specific color.

6.2 Flights Clustering and Visualization

The methodology proposed is currently tested by Airbus to assist in building health indicators for condition based maintenance. Health indicators are used for assessing the current state of some system and also for forecasting its future states and possible degradation (e.g. bleed, power systems, engine, APU, …). Airlines can be then informed that some systems should be maintained, so as to avoid any operational procedure at a given time horizon susceptible to cause e.g. delays, operational interruptions, etc …. The construction of a health indicator can be basically summarized as follows:

Collect health and usage data from various aircrafts (generally one has to consider similar ones). 2. 2.

Collect some operational events happening on these aircrafts due to some aircraft system errors (e.g. operational interruption, delays) 3. 3.

Identify anomalies in the health and usage data. 4. 4.

Identify some dependencies between health and usage data anomalies and operational events (by means of statistical hypothesis testing but also thanks to human expertise). 5. 5.

As soon as some dependencies are well identified, a health indicator is built.

The main barrier is the identification and the understanding of the anomalies. Different operational events are often recorded, corresponding to the degradation of different systems. Usually, a first stage of anomaly detection is performed, followed by a clustering of the anomalies listed for interpretation purpose. The major advantage of the approach proposed in this paper is that it directly provides a similarity measure between the anomalies. This strategy is illustrated by Fig. 3. The proposed method was applied on a dataset of $18553$ flights, each of which is characterized by $82$ parameters. In order to differentiate between anomalies corresponding to unusually large and small values, each feature is duplicated and each copy of a given feature is defined as the positive (resp. negative) value of the parameter above (resp. below) its mean value.

Fig. 3 and Fig. 4 display the clustering of $300$ ’extremal’ flights into $18$ groups, showing on the one hand the output of the spectral clustering applied to the similarity graph $w_{\hat{\bm{\theta}}}$ and on the other hand the underlying graph obtained with the same procedure as in Fig. 2.

6.3 A Real World Data Experiment with the Ground Truth

The shuttle dataset is available in the UCI repository, see Dheeru and Karra Taniskidou, (2017) (training and test datasets are merged here), $9$ numerical attributes and $7$ classes are observed. Class $1$ representing more than $80\%$ of the dataset, since our goal is to cluster rare and extreme events, instances from all classes but $1$ are analyzed, leading to a sample size equal to $12414$. The number of extreme points considered is denoted by $n_{0}$ here. We compare our approach to the $K$-means algorithm and the spectral clustering algorithm as implemented in Pedregosa et al., (2011). The number of clusters that we fix in advance to run each of these algorithms is denoted by $n_{cluster}$. The performance of each approach is evaluated by computing the purity score:

[TABLE]

where $n_{i,c}$ is the number of elements of class $c$ in the cluster $i$. As shown by Table 4, the purity score produced by the anomaly clustering technique promoted in this paper is always equal to or higher than those obtained by means of the other algorithms.

7 Conclusion

Because extreme values (viewed as anomalies here) cannot be summarized by simple meaningful summary statistics such as local means or modes/centroids, clustering and dimensionality reduction techniques for such abnormal observations must be of very different nature than those developed for analyzing data lying in high probability regions. This paper is a first attempt to design a methodology fully dedicated to the clustering and visualization of anomalies, by means of a statistical mixture model for multivariate extremes that can be interpreted as a noisy version of the angular measure, which distribution on the unit sphere exhaustively describes the limit dependence structure of the extremes. Mixture component are identified here with specific sub-simplices forming the support of the angular measure. Considering synthetic and real datasets, we also provide empirical evidence of the usefulness of (graph-based) techniques that can be straightforwardly implemented from the framework we developed.

Acknowledgements

This work was supported by a public grant as part of the Investissement d’avenir project, reference ANR-11-LABX-0056-LMH, LabEx LMH.