<figure><img src="image_1.png" /></figure>

HAL
open science

# Functional anomaly detection: a benchmark study

Guillaume Staerman, Eric Adjakossa, Pavlo Mozharovskyi, Vera Hofer, Jayant Sen Gupta, Stéphan Clémençon

<figure><img src="image_2.png" /></figure>

## To cite this version:

Guillaume Staerman, Eric Adjakossa, Pavlo Mozharovskyi, Vera Hofer, Jayant Sen Gupta, et al. Functional anomaly detection: a benchmark study. International Journal of Data Science and Analytics, 2023, 16 (1), pp.101-117. 10.1007/s41060-022-00366-5. http://dx.doi.org/10.1007/s41060-022-00366-5

HAL Id: hal-04268013

https://telecom-paris.hal.science/hal-04268013v1

Submitted on 1 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Functional Anomaly Detection: a Benchmark Study

Guillaume Staerman¹*, Eric Adjakossa¹, Pavlo Mozharovskyi¹, Vera Hofer², Jayant Sen Gupta³ and Stéphan Clémençon¹

¹*LTCI, Télécom Paris, Institut Polytechnique de Paris, France.

²Department of Operations and Information Systems, University of Graz, Austria.

³ Airbus AI Research, Toulouse, France.

*Corresponding author(s). E-mail(s): guillaume.staerman@telecom-paris.fr;

## Abstract

The increasing automation in many areas of the Industry expressly demands to design efficient machine-learning solutions for the detection of abnormal events. With the ubiquitous deployment of sensors monitoring nearly continuously the health of complex infrastructures, anomaly detection can now rely on measurements sampled at a very high frequency, providing a very rich representation of the phenomenon under surveillance. In order to exploit fully the information thus collected, the observations cannot be treated as multivariate data anymore and a functional analysis approach is required. It is the purpose of this paper to investigate the performance of recent techniques for anomaly detection in the functional setup on real datasets. After an overview of the state-of-the-art and a visual-descriptive study, a variety of anomaly detection methods are compared. While taxonomies of abnormalities (e.g. shape, location) in the functional setup are documented in the literature, assigning a specific type to the identified anomalies appears to be a challenging task. Thus, strengths and weaknesses of the existing approaches are benchmarked in view of these highlighted types in a simulation study. Anomaly detection methods are next evaluated on two datasets, related to the monitoring of helicopters in flight and to the spectrometry of construction materials namely. The benchmark analysis is concluded by a recommendation guidance for practitioners.

**Keywords:** Anomaly detection, Functional data, Depth statistics, Taxonomy of anomalies, Benchmark study

## 1 Introduction

The infatuation with machine learning is spreading to nearly all fields (e.g. science, transportation, energy, medicine, security, banking, insurance, commerce) as the ubiquity of sensors (e.g. IoT) make more and more data at disposal with an ever finer granularity. An abundance of new applications such as the health monitoring of complex infrastructures (e.g. aircrafts, energy networks) together with the availability of massive data samples (Big Data) and technological bricks for information acquisition and access (e.g. sensor

networks, IoT, distributed file systems) and computation (e.g. infrastructures for massively parallelized/distributed computation, on-line processing) has put pressure on the scientific community to develop new artificial intelligence methods and algorithms, following the example of successful applications such as computer vision, machine-listening or machine language translation. In particular, (unsupervised) anomaly detection remains a very active research topic. It can be stated in an informal manner as the identification of data
points $x$ that are suspicious due to their significant difference from the vast majority of observed data. The goal pursued is the same as in binary classification, to build a classification rule $g$ mapping the input space, $\mathcal{X}$ say, to $\{-1, +1\}$ namely: an observation $x \in \mathcal{X}$ is then considered as abnormal/anomalous, i.e. not distributed according to the same probability measure $P$ as that of the 'normal' (to be understood as 'not abnormal', with no reference to the Gaussian distribution here) observations when $g(x) = +1$, and as normal when $g(x) = -1$. However, in contrast to the supervised situation where a binary label is assigned to each training observation indicating whether it is normal, only unlabeled data, supposedly normal for the most part, are available for learning the decision rule $g$ in the unsupervised case. With the deployment of sensors monitoring their operation nearly continuously, machine-learning is expected to provide solutions for predictive maintenance of sophisticated systems, such as electricity grids or aircrafts, facilitating the early detection of "weak" signals that announce breakdowns, and serving to plan the replacement of components before their probable failure. Beyond the control of the false alarm rate, the challenge now essentially consists in fully exploiting the information collected, taking often the form of measurements of a physical variable sampled at a very high frequency. In this case, the observations cannot be treated as multivariate data and a functional approach is required. The (unsupervised) anomaly detection issue in the multivariate setting (when $\mathcal{X} \subset \mathbb{R}^d$ with $d \ge 1$) is well-documented in the literature and a large variety of dedicated techniques (way too numerous to be listed exhaustively here) have been proposed and investigated. The vast majority of the heuristics considered consist in turning supervised learning methods into unsupervised approaches, where label is replaced by rarity, anomalies being supposed to be rare by definition (e.g. SVM becoming one-class SVM, classification trees becoming isolation trees). In contrast, there has not been as much attention paid to the functional situation until now, see e.g. [1] or [2] and the references therein. Functional data are available in the form of functions, images and shapes or more general objects [3], their statistical analysis, referred to as functional data analysis (FDA in abbreviated form) has received much interest in

the last two decades, see e.g. [4], [5], [3]. In its simplest form, a functional dataset is a collection of curves observed on a grid $t_1 < \dots < t_p$, generally supposed to possess some smoothness properties to allow their recovery by means of appropriate interpolation/approximation schemes with a controlled error, see e.g. [6] for an account of FDA techniques. Though a very rich information can be carried by functional data, the downsides of FDA are challenges for designing feasible numerical procedures and for establishing a sound validity theoretical framework alike. Of course, anomaly detection for functional data can be reduced to the implementation of its counterpart in the multivariate case by means of straightforward dimensionality reduction techniques. Indeed, one may first project the observed curves onto a subspace of finite dimension by keeping the most informative components of a Kahrunen-Loève decomposition (i.e. those with largest empirical variance) or by truncating their expansion in a Hilbertian basis of reference, using a flexible dictionary of functions (e.g. Fourier basis, wavelets). This step is often referred to as the *filtering* stage. Next, any anomaly detection method designed for multivariate observations can be used for analyzing the (parsimonious) filtered data thus obtained, see, e.g. [5]. Though easy to implement, such a two-step strategy is unfortunately dramatically conditioned by the finite-dimensional representation method chosen, causing the loss of patterns possibly crucial to discriminate between normal observations and abnormal observations of certain types. Indeed, a taxonomy of functional outliers in [7] distinguish between *isolated outliers*, which manifest an outlying behavior during a very short time interval, like spikes or peaks in univariate curves, and *persistent outliers*, which exhibit a deviation from the normal behavior on a large part of the domain and consist of *shift outliers* (same shape as the majority of the observations, but moved away), *shape outliers* (curves with a shape which differs from that of the majority of the observations) or *amplitude outliers* (curves with a scale/range which differs from that of the majority of the observations). Hence, a certain finite-dimensional representation tailored to the detection of a specific type of anomaly, it may completely fail in detecting the other types. Because there exists a variety of functional anomalies, robust statistical methods such as those based
on depth statistics have recently received attention, see e.g. [8] or [9]. The depth of a point within a data cloud was first introduced in [10], in order to define a notion of multivariate median in $\mathbb{R}^d$ with $d \ge 2$. Afterwards, so as to measure the depth of an arbitrary point $x$, [11] considered hyperplanes through $x$ and determined its depth by the smallest portion of data that are separated by such a hyperplane. Since then, this idea has proved to be very fruitful and has lead to a rich statistical methodology, still in progress, in particular with the design of more general nonparametric depth statistics. In a nutshell, the depth of a data point describes its *outlyingness*, or *centrality* inversely, relatively to the data cloud, and thus defines a pre-order on the feature space: the smaller its depth, the likelier a point can be considered as an outlier. For a dataset at hand, in the majority of cases computation of depth of a point consists in optimising a specific function: e.g. minimization of the directional outlyingness, since an observation is considered abnormal if it is abnormal at least in one direction, while normal observations are those that are normal in all directions. Although depth functions can be computed for empirical distributions only from a practical perspective, the consistency results documented in the literature provide guarantee that statistical versions may provide accurate estimates of the depth of the true distribution of the data. Extension of these methods and results to data in functional spaces have been recently the subject of attention, see e.g. [12], with applications to functional anomaly detection in particular, refer to e.g., [7] or [13].

It is the goal of this article to investigate the performance of recent techniques for functional anomaly detection and compare their accuracy with that of simpler approaches, based on a preliminary dimensionality reduction, standing as natural competitors. In particular, specific attention is paid to those that are based on functional depth statistics or that extend multivariate methods by avoiding the filtering step. A benchmark study comparing the merits of the methods considered here regarding various metrics of reference is thus presented on aeronautics data gathered by Airbus and spectrometry measurements of sedimentary material collected by the Geological Survey of Austria for quality assessment on mining sites of Austria. Specifically, the aeronautics dataset consists of one-minute-sequences of

accelerometer data measured at a 1024Hz frequency. Airbus dataset is divided into two parts: the training set composed of 1677 curves with no available labels that may contain 'abnormal' observations and the validation/test set composed of 2511 time-series with 1794 'normal' curves. In contrast, the test data are labeled, in order to evaluate the performance of the anomaly detection rules learned in the training stage. These measurements were made on test helicopters at various locations, in various angles, on different flights. The learning framework is unsupervised: in the experiment, all accelerometer data series at disposal for training automatically a classifier to detect abnormal changes are considered as normal. The spectrometry of rocks data consists of one dataset of 2096 curves with 60000 measurements. It represents materials of two types whose labels are available, with limestone being the desirable (normal) rock type and the intrusive (abnormal) cellular dolomite. The task is thus, given the reflectance spectrum (with noise subtracted and normalized with respect to the reference spectrum) of multiple samples of mined examples, separate those abnormal.

As revealed by the experimental analysis carried out, recent (depth-based) functional anomaly detection techniques significantly outperform traditional methods. Additional experiments based on simulation data also provide empirical evidence that their flexibility permits to detect functional anomalies of different types, avoiding the limitations due to the exploitation of a specific finite-dimensional representation of the data.

The paper is organized as follows. Section 2 gives a general overview of recent anomaly detection methods for functional data, and discusses both their strengths and weaknesses. In Section 3, the performance metrics used for measuring the accuracy of anomaly detection rules learned in an unsupervised manner on (labeled) test data are described, the experiments on synthetic/real data and the results obtained are presented at length in Section 4. Finally, some concluding remarks are collected in Section 5.
# 2 Anomaly Detection for Functional Data

We start with recalling briefly the rationale behind classic approaches to functional anomaly detection and next describe more recent techniques dedicated to this task, coping directly with the functional nature of the observations.

## 2.1 Classic Approaches - Filtering Techniques

In the unsupervised learning framework, no label indicating whether a training observation is anomalous or not is available. Hence, anomalies should be identified in an automatic way by learning the 'normal' behavior, that of the vast majority of the observations, and considering those differing significantly from it as 'abnormal'. Logically, anomalies are rare in the data and thus fall in 'low density' regions: anomaly detection thus boils down to identifying the 'tail' of the distribution.

In the standard multivariate setup, i.e. when the observations can be represented as points in $\mathbb{R}^d$, a number of methods have been developed in the statistical literature. When the data generating process is known, a model based anomaly detection can be used of course. For example, if the distribution modelling the normal behavior is supposed to be Gaussian, observations outside a certain ellipsoid (corresponding to an isodensity curve) are considered as abnormal. In order to recover it, a robust mean and covariance estimation should be performed, see e.g. [14]. However, in most practical situations, the probability distribution of the normal observations is unknown and nonparametric approaches like Minimum Volume Set estimation (MV-set) [15, 16], One-Class Support Vector Machines (OCSVM) [17], Local Outlier Factor (LOF [18]) or Isolation Forest (IF) [19, 20] are widely used. Most of these algorithms produce, under appropriate conditions, estimates of upper-level sets of the density of the normal observations, tagging as anomalies the observations lying outside these regions. The concept of statistical depth, originally introduced by [10] (refer also to [21] for a list of desirable properties a depth function should fulfill), permits to extend the notion of quantile to distributions $P$ on $\mathbb{R}^d$ when $d \ge 2$ and can also be used for

the purpose of anomaly detection. By defining a center-outward ordering of points in the support of a probability distribution $P$, a depth function $D(\cdot | P) : \mathbb{R}^d \to \mathbb{R}^d$ serves to quantify the centrality of any point $x \in \mathbb{R}^d$: points with highest depth values are ideally those near the 'center' of the mass, while observations with depth values less than a certain critical threshold (i.e. outside a depth region $\{x \in \mathbb{R} : D(\cdot | P) \ge \tau\}$ with $\tau > 0$) are considered as abnormal, see e.g. [22].

<figure><img src="image_4.png" /><figcaption>Poorly reconstructed function</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 1** Two functional observations (black curves) from the aeronautics dataset and their reconstruction (blue curves) using 10 first principal components of FPCA: a normal observation, well reconstructed (top and an abnormal one, poorly reconstructed (bottom)).</figcaption></figure>

In order to extend the application of the techniques mentioned above to functional data, a filtering step can be used. Assuming that the variable observed takes its values in a subset of a Hilbert space, the data are first projected onto a finite dimensional subspace defined by truncating their expansion in an orthonormal basis and any multivariate anomaly detection technique can be next fed with the filtered data. The basis considered may be picked in a dictionary (e.g. Fourier,
wavelets, splines) or learned from the data as in Functional Principal Component Analysis (FPCA in abbreviated form). In the ideal case where one knows in advance which type of anomalies one attempts to detect, accuracy can be optimized by tuning the parameters of the filtering stage (basis, subspace which the data are projected onto) accordingly. However, in most situations encountered in practice, there are no guarantees about the nature of future anomalies: if the finite-dimensional representation produced by the filtering procedure can be appropriate to detect certain types of anomalies, it may completely fail in detecting the other types, see Figure 1. To remedy this drawback, certain methods tailored to functional data have been recently developed.

Table 1 presents the F1-Score, Average Precision, Area Under the Receiver Operating Characteristic and Sensitivity (detailed later in Section 3.1) for the data mentioned above through application of three very popular anomaly detection techniques (Isolation Forest, Local Outlier Factor, One-Class Support Vector Machine) after a projection on a finite-dimensional subspace defined by means of the Haar basis. Although the results are barely satisfactory at first glance, in view of the high complexity of the data, the greater performance of methods targeting directly the functional nature of the data will be demonstrated in the analysis carried out in Section 4.2. While a fine tuning of the employed dimension reduction technique may possibly improve the result, this is far from simple in practice, insofar as it could require to implement an extremely large number of methods with many parameter configurations.

## 2.2 Coping with the Functional Nature of the Data

We now review recent functional depth based techniques, which may offer attractive practical alternatives to the traditional techniques previously recalled. Suppose that the random variable of interest $X$ modelling the normal behavior of the system monitored takes its values in a Hilbert space, say the space $L^2([0, 1])$ of square integrable functions on $[0, 1]$ for simplicity. Direct extension of multivariate data depth methods to this functional setting turns to be impractical because the resulting depth functions then vanish everywhere

<table><caption>**Table 1** Finite-dimensional algorithms on Haar-basis projection coefficients considered in performance comparison with F1-Score, Average Precision (AP), AUC and Sensitivity ($p_c$) for both datasets.</caption><thead><tr><th>Method</th><th>F1-Score</th><th>AP</th><th>AUC</th><th>p<sub>c</sub></th></tr></thead><tbody><tr><td colspan="5">Aeronautics dataset</td></tr><tr><td>IF</td><td>0.71</td><td>0.69</td><td>0.49</td><td>0.27</td></tr><tr><td>LOF</td><td>0.68</td><td>0.73</td><td>0.5</td><td>0.21</td></tr><tr><td>OCSVM</td><td>0.78</td><td>0.77</td><td>0.63</td><td>0.45</td></tr><tr><td colspan="5">Rocks dataset</td></tr><tr><td>IF</td><td>0.971</td><td>0.978</td><td>0.528</td><td>0</td></tr><tr><td>LOF</td><td>0.973</td><td>0.965</td><td>0.495</td><td>0.05</td></tr><tr><td>OCSVM</td><td>0.971</td><td>0.965</td><td>0.686</td><td>0</td></tr></tbody></table>

during the optimization, as pointed out in e.g. [23], due to the richness of the feature space. Further, a set of desirable properties is imposed on the defined functional depth statistic, which guarantees its usefulness in applications, see [24, 25] for their overview. To adjust for these requirements, various strategies have been proposed. In particular, the search of the depth function is restricted to features from a dictionary of interest only in [26], while integrating univariate depth functions over time is suggested in e.g. [27] or [7]. This second approach is often preferred in practice, due to its simplicity and (often substantially) lower computational burden.

**Integrated functional depths.** Let $D^{(1)}(\cdot | \cdot)$ be a univariate depth measure, refer to e.g. [21] and [22] for an account of depth statistics. A typical example is the Tukey depth, see [10] and also [28] or [27]. For an i.i.d. sample $D_n = \{X_1, \dots, X_n\}$ of real valued observations, it is defined as:

$$D_{T}^{(1)}(x | D_{n}) = \min\{\hat{F}_{n}(x), 1 - \hat{F}_{n}(x^{-})\}, \quad (1)$$

where the càd-làg function $\hat{F}_n(x) = (1/n) \sum_{i=1}^{n} \mathbb{I}\{X_i \le x\}$ is the empirical cdf based on the dataset $D_n$, denoting by $\mathbb{I}\{\mathcal{E}\}$ the indicator function of any event $\mathcal{E}$. Now, given a sample of curves $C_n = \{X_1, \dots, X_n\}$ composed of independent copies of the functional observation $X(t)$ with $t \in [0, 1]$ valued in $L_2([0, 1])$, the integrated functional data depth function defined below provides a centrality measure on the infinite dimensional space $L^2([0, 1]):$
$$D(\mathbf{X} | \mathcal{C}_n) = \int_0^1 D^{(1)}(\mathbf{X}(t) | \mathcal{C}_n(t)) dt, \quad (2)$$

where $\mathcal{C}_n(t) = \{\mathbf{X}_1(t), \dots, \mathbf{X}_n(t)\}$ is the stamp of the functional sample at time point $t$. While being purely data-driven and robust by construction, the univariate data depth (1) is non-continuous, and equals 0 outside the range of the $X_i$'s. Another option is thus to consider the integrated functional depth based on the projection depth introduced in [21], which is given by:

$$D_P^{(1)}(x | \mathcal{D}_n) = \left( 1 + \frac{|x - \text{med}(X)|}{\text{MAD}(X)} \right)^{-1},$$

where med($X$) and MAD($X$) are the median and the median absolute deviation based on the univariate sample $\mathcal{D}_n$. While being positive everywhere, projection depth is symmetric around the median. More generally, the asymmetric version proposed by [7] is given by

$$D_{AP}^{(1)}(x | \mathcal{D}_n) = \frac{1}{1 + \text{AO}(x | \mathcal{D}_n)},$$

where the asymmetric outlyingness AO is defined as:

$$\text{AO}(x | \mathcal{D}_n) = \begin{cases} \frac{x - \text{med}(X)}{w_2(X) - \text{med}(X)} & \text{if } x > \text{med}(X), \\ \frac{\text{med}(X) - x}{\text{med}(X) - w_1(X)} & \text{if } x < \text{med}(X) \end{cases}$$

with

$$w_1(X) = \hat{q}_X(0.25) - \text{IQR}(X)1.5e^{-4\text{MC}(X)}$$

$$w_2(X) = \hat{q}_X(0.75) + \text{IQR}(X)1.5e^{3\text{MC}(X)},$$

where $\text{IQR}(X) = \hat{q}_X(0.75) - \hat{q}_X(0.25)$ is the interquantile range with $\hat{q}_X$ being the empirical quantile based on $\mathcal{D}_n$ and $\text{MC}(X)$ is the robust measure of skewness proposed by [29] defined as follows:

$$\text{med} \left( \frac{(X_j + X_i - 2 \text{med}(X))}{x_j - x_i} : (i, j) \in \mathcal{I} \right),$$

where $\mathcal{I} = \{(i, j) : i \neq j, X_i \leq \text{med}(X) \leq X_j\}$. In the subsequent empirical analysis, three integrated functional data depth functions are considered, calculated as the average Tukey depth (fT),

projection depth (fSDO) and asymmetric projection depth (fAO) over the definition domain, respectively.

<figure><img src="image_6.png" /><figcaption>Fig. 2 Illustration of the influence of the convex hull of a single function in the ACH depth (left, top and bottom, abnormal (red) observation increases the volume) and separation in the scalar product in FIF (right, top and bottom, abnormal (red) observation lies far from the normal data).</figcaption></figure>

## The Area of the Convex Hull (ACH) depth.

The area of the convex hull (ACH) depth, recently introduced by [9], is also considered in the present study. We point out that this functional depth function does not belong to the family of integrated depths (2) and exhibits high sensitivity beyond the convex envelope of the data. In short, the ACH depth quantifies the contribution of a given curve to the spread of the convex hull of a random (sub)sample of curves, see Figure 2 left. For a chosen degree $1 \le J \le n$, the ACH depth, denoted by $D_{ACH}(\mathbf{X} | \mathcal{C}_n)$ here, is defined as follows:

$$\frac{1}{\binom{n}{J}} \sum_{1 \le i_1 \le \dots \le i_j \le n} \frac{\lambda(\text{conv}(\text{graph}(\{\mathbf{X}_{i_1}, \dots, \mathbf{X}_{i_j}\})))}{\lambda(\text{conv}(\text{graph}(\{\mathbf{X}_{i_1}, \dots, \mathbf{X}_{i_j}\} \cup \{\mathbf{X}\})))},$$

where $\lambda$ means the Lebesgue measure on $\mathbb{R}^2$, $\text{conv}(A)$ is the convex hull of any subset $A$ of the plane and $\text{graph}(\{\mathbf{X}_1, \dots, \mathbf{X}_n\}) = \bigcup_{i=1}^n \text{graph}(\mathbf{X}_i)$ with the notation $\text{graph}(\mathbf{X}_i) = \{(t, y) : y = \mathbf{X}_i(t), t \in [0, 1]\}$. In [9] empirical evidence that even small values of $J$ ($J = 2$ for instance) are sufficient to define a robust
and sensitive functional depth that is relevant for anomaly detection, the parameter $J$, when picked in this range, having little impact on the pre-order induced by the resulting depth function. One may refer to [9] for a detailed description of this approach.

**Functional Isolation Forest.** Recently, [2] extended the popular Isolation Forest approach, see [19], to the functional framework. A Functional Isolation Forest (FIF) is formed by a collection of functional isolation trees (F-ittrees), constructed each (by a sequence of random splits) from a subsample (of size $\psi$) of $C_n$. The abnormality score of the observation $\mathbf{X}$ is then computed as a monotone decreasing transformation of the average depth of $\mathbf{X}$ over the trees. The main idea here is that, since the splits are purely random, a very differing observation will be cut-off (isolated) from $C_n$ with higher probability (and thus on less deep levels of the F-ittrees) than those similar to the majority of the observed curves. The construction of the F-ittrees is based on a pre-determined dictionary $\mathbf{D}$ that can contain both deterministic and/or stochastic functions that capture relevant properties of the data, which can also be a subset of $C_n$. Before each random univariate split, all observations of the F-itlee are projected on a line spanned by a random element of $\mathcal{D}$, see Figure 2 right. The choice of a suited dictionary thus plays a crucial role in the construction of the FIF score. This projection is calculated by means of the scalar product designed to account for both location and shape anomalies, i.e. for any $\alpha \in [0, 1]$:

$$ \langle \mathbf{X}, d \rangle := \alpha \times \frac{\langle \mathbf{X}, d \rangle_{L_2}}{\|\mathbf{X}\|\|d\|} + (1 - \alpha) \times \frac{\langle \mathbf{X}', d' \rangle_{L_2}}{\|\mathbf{X}'\|\|d'\|}, $$

where $d \in \mathbf{D}$, $X'$ and $\|X\|$ are respectively the first derivative and the $L_2$-norm of $X$ in $L_2([0, 1])$.

## 3 A Preparatory Simulation Study

This section is devoted to empirical analysis of the performance of the functional anomaly detection techniques whose mechanics have been briefly described in Section 2. Their accuracy is investigated from simulated data inspired from a real

dataset collected by Airbus, composed of one-minute sequences of accelerometer data measured on helicopters. As a first go, we recall the standard performance metrics commonly used.

### 3.1 Performance Metrics in Anomaly Detection

Although the learning procedure does not rely on data of the same nature, the anomaly detection problem can be formulated in the same probabilistic framework as binary classification, the flagship problem in statistical learning. In the standard setup, the binary random variable $Y$ indicates the occurrence of an anomaly: the label is positive, i.e. $Y = +1$, when an anomaly occurs, and negative, i.e. $Y = -1$, otherwise. The random variable $\mathbf{X}$, taking its values in the feature space $\mathcal{X}$, models the measurement at disposal to predict $Y$. The goal pursued is to build an anomaly scoring function $s: \mathcal{X} \to \mathbb{R} \cup \{+\infty\}$ in an unsupervised manner (without observing $Y$), so that, ideally, the larger $s(\mathbf{X})$, the likelier an anomaly occurs, i.e. the more probable is $Y = +1$. A decision to raise an alarm can be then built by thresholding the scoring function $s(\mathbf{x})$ at a critical level, ruling the trade-off between errors of type I and type II. Equipped with this notation, decreasing transforms of a depth function w.r.t. $\mathbf{X}$'s marginal distribution provide anomaly scoring functions in a natural fashion. Precisely, when using data depths in the following study, the transformation $1 - D(.)$ is performed to rescale them as an anomaly score.

**ROC analysis.** The golden standard to quantify theoretically the accuracy of an anomaly scoring function $s$ is the PP-plot of the false positive rate $\text{vs}$ the true positive rate, namely $t \in \mathbb{R} \mapsto (\mathbb{P}\{s(\mathbf{X}) \ge t | Y = -1\}, \mathbb{P}\{s(\mathbf{X}) \ge t | Y = +1\})$, referred to as the ROC curve (standing for Receiver Operator Characteristic curve), see e.g. [30]. The higher the curve, the more accurate the anomaly scoring function. A simple Neyman-Pearson argument shows that optimal scoring functions are increasing transforms of the likelihood ratio $\Psi(\mathbf{X}) := (dF_+/dF_-)(\mathbf{X})$, denoting by $F_\sigma$ the conditional distribution of $\mathbf{X}$ given $Y = \sigma 1, \sigma \in \{-, +\}$: their ROC curve dominating everywhere the ROC curve of any other anomaly scoring function. For this reason, this functional performance measure is generally summarized by
the Area Under the ROC curve (AUC in abbreviated form), a popular scalar criterion that can be classically interpreted as the rate of concordance of pairs: AUC(s) = P{s(X) > s(X') | Y = +1, Y' = -1} + P{s(X) = s(X') | Y = +1, Y' = -1}/2, where (X', Y') denotes an independent copy of the pair (X, Y).

**PR analysis.** Alternatively, one may evaluate the accuracy of any score by plotting the precision-recall (shortly PR) curve, namely $t \mapsto (P\{s(X) \ge t | Y = +1\}, P\{Y = +1 | s(X) \ge t\})$. The higher its PR curve, the more accurate an anomaly scoring function. Of course, as may be immediately shown by means of the Bayes formula, PR and ROC curves are in one to one correspondence, see e.g. [31]. Like the ROC curve, the PR curve may be summarized by the area under it, referred to as the Average Precision (AP).

When labeled data are available, it is possible to compute the performance measures recalled above, replacing the probabilities involved by their statistical counterparts, in order to assess the accuracy of any anomaly scoring function candidate $s(x)$. However, in unsupervised anomaly detection, the scoring function $s(x)$ cannot be learned using labeled training data, in contrast to classification or bipartite ranking: only 'negative' observations, i.e. an i.i.d. sample drawn from (a possibly noisy version of) distribution $F_-$, are available in the training stage. Hence, the learning task cannot be achieved by optimizing empirical versions of the aforementioned criteria, which makes it extremely challenging.

## 3.2 Simulating Anomalies of Specific Types

Datasets containing various types of anomalies are usually more challenging to analyze. As a first go, we start by investigating to which extent the techniques recalled above may permit to detect simulated anomalies of well-identified types according to the usual taxonomy, see, e.g., [7]. Figure 3 illustrates the four types of anomalies addressed in detail in these experiments: isolated, magnitude (of two different kinds) and shape anomalies.

To reproduce a controlled version of each type of anomalies, four datasets have been built from a collection of 1794 'normal' functional observations from the validation dataset collected by Airbus. Each functional observation corresponds

to accelerometer data measured on helicopters at a 1024 Hz frequency over time windows of 1 minute: the curves $X = (X(t))_{t \in [0,1]}$ are built by means of an affine interpolation of the 61440 sampled points. One per anomaly type, four datasets have been constructed by adding a specific contamination to 5% of these 'normal' observations, drawn uniformly at random. Cases when 1%, 2% 3% and 4% are added in the Appendix B for completeness and show similar behavior of methods than for 5%. The four contamination models defined below, are used to generate independent curves $Y$ (independently from the original dataset) that are next added to the selected above 1794 'normal' observations $X$. By $U([a, b])$ is meant the uniform distribution on the interval $[a, b]$, while $\delta_u$ denotes the Dirac mass at point $u$.

**Model 1 (Isolated Anomalies)**

$Y(t) = \varepsilon u_1 I\{t = \tau\}$, where $u_1 \sim U([3, 4])$ and $\varepsilon \sim (1/2)(\delta_{-1} + \delta_1)$ are independent random variables, with $\tau$ being the time at which the isolated anomaly occurs that is chosen randomly in a uniform manner among the set of sampling points, independently of $u_1$ and $\varepsilon$.

**Model 2 (Magnitude Anomalies I)**

$Y(t) = u_2$, with $u_2 \sim U([-12, -15])$.

**Model 3 (Magnitude Anomalies II)**

$Y(t) = u_3 I(t \in I)$, where $u_3 \sim U([0, 15])$ and $I$ is a subinterval of $[0, 1]$ of length 1/10 whose location is chosen uniformly at random, independently of $u_3$.

**Model 4 (Shape Anomalies)**

$Y(t) = \sin(2\pi u_4 t)$, where $u_4 \sim U([0.2, 2])$.

Simulated anomalies, together with a small subset of 'normal' data, are illustrated in Figure 4.

<table><caption>**Table 2** Methods considered in performance comparison with the sensitivity ($p_c$) and the Area Under the Receiver Operating Characteristic (AUROC) for the four simulated models.</caption><thead><tr><th>Anomaly type</th><th></th><th>IF</th><th>LOF</th><th>OCSVM</th></tr></thead><tbody><tr><td rowspan="2">Isolated</td><td>p<sub>c</sub></td><td>0</td><td>0</td><td>0</td></tr><tr><td>AUC</td><td>0.41</td><td>0.25</td><td>0.44</td></tr><tr><td rowspan="2">Magnitude I</td><td>p<sub>c</sub></td><td>1</td><td>0.48</td><td>1</td></tr><tr><td>AUROC</td><td>1</td><td>0.97</td><td>1</td></tr><tr><td rowspan="2">Magnitude II</td><td>p<sub>c</sub></td><td>0</td><td>0</td><td>0</td></tr><tr><td>AUC</td><td>0.54</td><td>0.03</td><td>0.7</td></tr><tr><td rowspan="2">Shape</td><td>p<sub>c</sub></td><td>0</td><td>0.48</td><td>0</td></tr><tr><td>AUC</td><td>0.67</td><td>0.97</td><td>0.67</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>**Fig. 3** Examples for each of the identified types of anomalies in the aeronautics dataset. In order from top to bottom and from left to right: each time one isolated, magnitude, magnitude/shape and shape anomaly, with a subsample of normal data. Grey curves are normal data while red curves are anomalies.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 4** Examples for each of four types of simulated anomalies. In order from top to bottom and from left to right: each time one isolated, magnitude, magnitude/shape and shape anomaly, with a subsample of normal curves (out of 1794). Grey curves are normal data while red curves are simulated anomalies.</figcaption></figure>
<table><caption>**Table 3** Methods considered in performance comparison with the sensitivity ($p_c$) and the Area Under the Receiver Operating Characteristic (AUROC) for the four simulated models with 5% of added anomalies.</caption><thead><tr><th rowspan="2">Methods</th><th colspan="2">Isolated</th><th colspan="2">Magnitude I</th><th colspan="2">Magnitude II</th><th colspan="2">Shape</th></tr><tr><th>p<sub>c</sub></th><th>AUC</th><th>p<sub>c</sub></th><th>AUC</th><th>p<sub>c</sub></th><th>AUC</th><th>p<sub>c</sub></th><th>AUC</th></tr></thead><tbody><tr><td>FIF</td><td>0</td><td>0.21</td><td>0.93</td><td>0.99</td><td>0</td><td>0.30</td><td><strong>0.63</strong></td><td><strong>0.98</strong></td></tr><tr><td>fAO</td><td>0</td><td>0.44</td><td><strong>1</strong></td><td><strong>1</strong></td><td>0</td><td>0.54</td><td>0</td><td>0.66</td></tr><tr><td>fbd</td><td>0</td><td>0.44</td><td><strong>1</strong></td><td><strong>1</strong></td><td>0</td><td>0.54</td><td>0</td><td>0.68</td></tr><tr><td>fSDO</td><td>0</td><td>0.42</td><td><strong>1</strong></td><td><strong>1</strong></td><td>0</td><td>0.43</td><td>0</td><td>0.77</td></tr><tr><td>fT</td><td>0</td><td>0.43</td><td><strong>1</strong></td><td><strong>1</strong></td><td>0</td><td>0.44</td><td>0</td><td>0.71</td></tr><tr><td>ACH</td><td>0</td><td><strong>0.63</strong></td><td>0.48</td><td>0.97</td><td><strong>0.80</strong></td><td><strong>0.99</strong></td><td>0</td><td>0.56</td></tr><tr><td>Outliergram</td><td>0</td><td>0.55</td><td><strong>1</strong></td><td><strong>1</strong></td><td>0</td><td>0.54</td><td>0</td><td>0.47</td></tr><tr><td>MS + IF</td><td>0</td><td>0.05</td><td>0.66</td><td>0.98</td><td>0</td><td>0.70</td><td>0.33</td><td>0.74</td></tr><tr><td>FOM (fSDO) + IF</td><td>0</td><td>0</td><td>0.85</td><td>0.99</td><td>0.64</td><td>0.99</td><td>0.06</td><td>0.80</td></tr><tr><td>FOM (fAO) + IF</td><td>0</td><td>0.14</td><td>0.81</td><td>0.99</td><td>0.55</td><td>0.98</td><td>0.02</td><td>0.87</td></tr><tr><td>FPCA + IF</td><td>0</td><td>0.11</td><td>0</td><td>0.91</td><td>0</td><td>0.71</td><td>0.45</td><td>0.97</td></tr><tr><td>FPCA + LOF</td><td>0</td><td>0.5</td><td>0</td><td>0.16</td><td>0</td><td>0.38</td><td>0</td><td>0.79</td></tr><tr><td>FPCA + OC</td><td>0</td><td>0.04</td><td>0</td><td>0.93</td><td>0</td><td>0.77</td><td>0.3</td><td>0.96</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Fig. 5** ROC curves of the four best methods for the simulation study.</figcaption></figure>
## 3.3 Naive Approaches - Sampled Curves Viewed as Multivariate Data

Being the most straightforward idea, and still frequently employed in functional anomaly detection, direct consideration of (discretized) functional data in the multivariate space (i.e., $\mathbb{R}^p$) of their measurements can be seen as a first naïve approach. It is important to notice that the sampling design should not vary with the curve/signal and preferably correspond to regularly spaced points in the observation domain, so that no dimensions are disadvantaged. If not, one can resort to importance-weighting techniques though, e.g., giving more weights to coordinates with higher marginal variance as it is suggested, e.g., in [27]. When implementing such a naïve strategy to deal with functional observations, specific attention must be paid to the possibly very large dimension of the data (compared to the size of the population sample), due to the high frequency character of the measurements, as it is the case for the aeronautics dataset considered in this paper where $61440 = p \gg n = 1677$. We have applied three of the most widely used methods for multivariate anomaly detection, namely Isolation Forest (IF), Local Outlier Factor (LOF), and One-Class Support Vector Machine (OCSVM) to the four settings described in Section 3.2.

Their sensitivity and Area Under Receiver Operating Characteristic (AUC) are reported in Table 2.

## 3.4 Results and Discussion

We now consider the set of methods which are in the center of attention of the current work. These treat the data directly in their original functional space, and—as we shall see right below—prove beneficial for anomaly detection.

The performance of these unsupervised methods is evaluated on the four simulated datasets described in Section 3.2 using the sensitivity (portion of correctly identified anomalies, $p_c$) and AUC as metrics. All parameters of the used algorithms are set to their default values (as it is pre-defined in the corresponding software packages). Thus, fSDO, fAO, fbd, FOM are available in the mrfDepth R-package [32]; the outliergram is available in the roahd R-package [33]; IF, LOF and

OC are available in `sklearn python` library [34]; Magnitude-Shape plot (MS) is available in the `scikit-fda python` library <fnref n="1" />; Functional Isolation Forest (FIF) and ACH open-source python codes are available under the following link <fnref n="2" />; fT can be easily coded from scratch.

Results of the simulation study are displayed in Table 3. For completeness, the ROC curves of the four best methods for each contamination setting are displayed in Figure 5. As expected, the score drastically varies across contamination models and anomaly detection methods. Isolated anomalies (especially short ones) of the Contamination Model 1 are difficult to detect with projections on most bases as well as by integrating depths, whilst ACH is sensitive to this kind of anomalies. Magnitude (especially type I) anomalies are known to be easier to detect and a number of methods (fAO, fbd, fSDO, fT and outliergram) perform well by managing to detect all of them. The difficulty that differentiates Magnitude II anomalies (from those in Magnitude I) is that the anomalies are expressed only for a subset of time points. This impedes many methods from detecting this kind of anomalies, while ACH seems to perform best among differing results, most probably due to slight resemblance of Magnitude II anomalies with the isolated ones. Shape anomalies is the least identifiable type, and FIF delivers better performance than other methods while taking into account both location and slope of the functional curves (due to the employed Sobolev-type metric).

## 4 Benchmarking Methods for Functional Anomaly Detection using Real Data

This section is devoted to the empirical analysis of the performance of the functional anomaly detection techniques whose mechanics have been briefly described in Section 2. Their accuracy, previously investigated based on artificially contaminated data, shall be now benchmarked using real labeled datasets.

<footnote n="1">¹https://github.com/GAA-UAM/scikit-fda</footnote>

<footnote n="2">²https://github.com/GuillaumeStaermanML</footnote>
<figure><img src="image_10.png" /><figcaption>**Fig. 6** Four visualization plots from Airbus data. In order from top to bottom and from left to right: FPCA, Magnitude-Shape plot with Integrated Tukey depth, Functional Outlier map with fAO and Functional Outlier map with fSDO. Black points corresponds to normal curves, red points to anomalies.</figcaption></figure>

<table><caption>**Table 4** Methods considered in performance comparison with F1-Score, Average Precision (AP), AUC and sensitivity ($p_c$) for Airbus and Rocks data. Bold numbers correspond to the best result.</caption><thead><tr><th rowspan="2">Method</th><th colspan="4">Airbus data</th><th colspan="4">Rocks data</th></tr><tr><th>F1-Score</th><th>AP</th><th>AUC</th><th>p<sub>c</sub></th><th>F1-Score</th><th>AP</th><th>AUC</th><th>p<sub>c</sub></th></tr></thead><tbody><tr><td>FIF</td><td><strong>0.81</strong></td><td><strong>0.88</strong></td><td><strong>0.76</strong></td><td>0.52</td><td>0.974</td><td><strong>0.992</strong></td><td>0.772</td><td>0.10</td></tr><tr><td>fAO</td><td>0.78</td><td>0.77</td><td>0.63</td><td>0.45</td><td><strong>0.989</strong></td><td>0.991</td><td><strong>0.833</strong></td><td><strong>0.58</strong></td></tr><tr><td>fbd</td><td>0.78</td><td>0.77</td><td>0.63</td><td>0.45</td><td>0.977</td><td>0.988</td><td>0.751</td><td>0.19</td></tr><tr><td>fSDO</td><td>0.78</td><td>0.77</td><td>0.63</td><td>0.45</td><td>0.972</td><td>0.980</td><td>0.555</td><td>0</td></tr><tr><td>fΓ</td><td>0.78</td><td>0.77</td><td>0.63</td><td>0.45</td><td>0.984</td><td>0.989</td><td>0.780</td><td>0.43</td></tr><tr><td>ACH</td><td>0.77</td><td>0.77</td><td>0.62</td><td>0.44</td><td>0.972</td><td>0.961</td><td>0.280</td><td>0</td></tr><tr><td>Outliergram</td><td>0.71</td><td>0.76</td><td>0.55</td><td>0.28</td><td>0.974</td><td>0.981</td><td>0.66</td><td>0.03</td></tr><tr><td>MS + IF</td><td>0.80</td><td>0.76</td><td>0.64</td><td>0.51</td><td>0.972</td><td>0.981</td><td>0.601</td><td>0</td></tr><tr><td>FOM (fSDO) + IF</td><td><strong>0.81</strong></td><td>0.76</td><td>0.66</td><td>0.53</td><td>0.972</td><td>0.978</td><td>0.530</td><td>0.02</td></tr><tr><td>FOM (fAO) + IF</td><td>0.80</td><td>0.77</td><td>0.65</td><td>0.51</td><td>0.984</td><td>0.991</td><td>0.804</td><td>0.448</td></tr><tr><td>FPCA + IF</td><td><strong>0.81</strong></td><td>0.80</td><td>0.70</td><td>0.53</td><td>0.972</td><td>0.971</td><td>0.446</td><td>0</td></tr><tr><td>FPCA + LOF</td><td>0.72</td><td>0.73</td><td>0.52</td><td>0.31</td><td>0.972</td><td>0.969</td><td>0.445</td><td>0.02</td></tr><tr><td>FPCA + OC</td><td><strong>0.81</strong></td><td>0.79</td><td>0.70</td><td><strong>0.54</strong></td><td>0.972</td><td>0.971</td><td>0.463</td><td>0</td></tr></tbody></table>
We analyze the complete Airbus dataset (i.e. including normal and abnormal accelerometer data) as well as data arising from the spectroscopy of rocks with wavelengths of light source ranging from 382 to 930 nanometers corresponding to visible-IR spectrum.

The complete Airbus dataset is split into an unlabeled training dataset (1677 sampled curves observed at 61440 equally spaced time points) used for the learning purposes and a labeled test dataset (2511 observations among which 717 are labeled as abnormal) to evaluate performance using the metrics described in Section 3.1. The second dataset contains linearly interpolated spectroscopy measurements of construction materials (mined rocks) from different locations and mining sites. Precisely, it includes 2038 limestones and 58 cell limes of 600 wavelength measures. Since the measured spectra are very noisy, and with the task of the current work being rather comparative (than absolute) performance of the methods, we modified the original dataset by removing most difficult to detect anomalies in a supervised manner (which only underlines the difficulty of the real-data anomaly detection task). More precisely, we removed 860 normal data that have (on an average) same values as the average value of anomalies and that can not be distinguished by any of the algorithms. Such removal shall not prioritize any of the methods. We further perform a linear interpolation in order to obtain data with equidistant measurement points.

## 4.1 Visualization

While a simple plot of the rocks data provides relevant information on the data shape (see the Supplementary Material), the richness of the aeronautics curves makes it impossible. To get a first insight into the structure of the aeronautics dataset under study, we start with meaningful visualization. Recently, many graphical tools have been proposed in the literature, such as general-purpose functional highest density region plots [35], functional boxplots [36] or amplitude and phase boxplot displays [37], as well as those especially designed for anomaly detection such as the Outliergram [38], the Functional Outlier Map (FOM) [7, 39] or the Magnitude-Shape plot (MS) [40]. This last group, together with the

generic Functional Principal Component Analysis (FPCA) constitute our interest.

The Magnitude-Shape plot is two-dimensional outlyingness (or data depth)-based graphical tool, based on the outlyingness mean and variance—over the time domain—of the functional observation. For a sample of curves $C_n = \{X_1, ..., X_n\}$ observed on $\{t_1, ..., t_p\}$, the MS plot is then the scatter plot of the points $(MO_i, VO_i)_{1 \le i \le n}$ with the Mean Outlyingness (MO)

$$MO_i = \frac{1}{p} \sum_{j=1}^{p} O(X_i(t_j) | C_n(t_j)),$$

and the Variational Outlyingness (VO)

$$VO_i = \frac{1}{p} \sum_{j=1}^{p} (O(X(t_j) | C_n(t_j)) - MO_i)^2,$$

where the outlyingness itself is usually a monotone decreasing transform of data depth:

$$O(X_i(t_j), C_n(t_j)) = \frac{1}{D(X_i(t_j) | C_n(t_j)) - 1}.$$

The Functional Outlier Map (FOM) is similar to the MS plot in the case of univariate functional data, with the only difference in measuring relative instead of absolute variability adapting thus to the actual variability of the function, see [39]. It is defined as a plot of points

$$(MO_i, \frac{\sqrt{VO_i}}{1 + MO_i})_{1 \le i \le n}.$$

For explainability purposes, vizualisation tools are computed on the test set where labels are available. Visualization plots, displayed in Figure 6, allow to identify certain anomalies (red points correspond to labeled anomalies), while others substantially overlap with the normal data (black points). They reveal the difficulty of marking the entity of abnormal observations for the aeronautics data, and explain variations in the performance of different methods used in Section 4.2.
## 4.2 Benchmark Study on
Aeronautics and Rocks Data

In Section 3.2 only four types of anomalies were identified and studied in more detail, while many others remain that are not easy to associate with any existing taxonomy, which is often the challenge of real-world data. While many methods are designed for aiming at certain types of anomalies, clearly no universality in detecting abnormal observations of different types can be generally expected.

To account for multiplicity of possible goals, we use several performance metrics: F1-Score, Average Precision (AP), AUC, and sensitivity $p_c$; see Section 3.1 for more details.

Table 4 displays the results. Methods perform differently and there is no general winner. First observation is the evidence of complexity of the real-world aeronautics dataset, which leads to non-perfect results across all the considered methods. This is also due to the fact that test data contains types of anomalies not present in the train data. Second, one should note that FIF appears to have very good performance, being best method in this benchmark in general. Even if its AUC of 76% is not high, it can be still seen as satisfactory provided existence of 29% of anomalies in the test data. Third, the depth-based methods indicate very stable results (mostly relatively satisfactory, except $p_c$), which comes from the fact that ordering of observations may be very similar for univarite functional curves due to (almost) coincidence of orderings for different univariate depths. While the rocks dataset was artificially simplified, the comparison remains similar and thus strengthens the conclusions.

Additionally, we display the ROC curves of the four best methods for both datasets (see Figure 7). The true positive rate is maximized by the FIF algorithm when the false positives are very close to zero which is crucial in many applications such as predictive maintenance in the case of airbus. This additional analysis makes FIF preferable to fAO in practice.

# 5 Conclusion and Perspectives

Due to its unsupervised nature, the anomaly detection problem is extremely challenging. Because the learning algorithms used to build a

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 7 ROC curves of the four best methods for the aeronautics (top) and rocks (bottom) data.</figcaption></figure>

decision rule in this context rely on 'normal' data solely, they mostly consist in identifying regions of the feature space where normal observations are less likely to fall in and cannot be based on empirical estimates of the performance metrics that shall be used afterwards to evaluate them when labels become available. This is even much more challenging in the case of observations valued in a vast functional space. In absence of prior knowledge about the types of the functional anomalies to be detected in the future, it is key to implement very flexible algorithms, exploiting the statistical information at disposal as far as possible. In this paper, we investigated the performance of recent anomaly detection techniques tailored to the functional framework such as Functional Isolation Forest and ACH depth, by means of real datasets, composed of sequences of accelerometer data measured on helicopters and spectrometry of construction materials. As confirmed by additional simulation studies, the benchmark revealed a clear advantage of these techniques compared to more traditional methods, relying on a preprocessing of the data (e.g. FPCA). These very encouraging results call for further applications
and extensions, to anomaly detection based on
multivariate time-series in particular, the behavior
of complex infrastructure being often continuously
monitored by several sensors and not just one.

Funding

This work has been funded by BPI France
in the context of the PSPC Project Expresso
(2017-2021). This project also received finan-
cial support from the initiative “Forschungspart-
nerschaften Mineralrohstoffe - ein strategischer
Forschungsschwerpunkt der Geologischen Bundes-
sanstalt”. The spectroscopic data of sedimentary
material was provided by the Geological Survey of
Austria.

References

[1] Peter J Rousseeuw and Mia Hubert. Anomaly detection by robust statistics. *Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery*, 8(2):e1236, 2018.

[2] Guillaume Staerman, Pavlo Mozharovskyi, Stéphan Clémençon, and Florence d'Alché Buc. Functional isolation forest. In *Proceedings of The 11th Asian Conference on Machine Learning*, pages 332–347, 2019.

[3] Jane-Ling Wang, Jeng-Min Chiou, and Hans-
Georg Müller. Functional data analysis.
*Annual Review of Statistics and Its Applica-
tion*, 3:257–295, 2016.

[4] James O. Ramsay and Bernard W. Silverman. *Functional Data Analysis*. Springer-Verlag, New York, 2005.

[5] Frédéric Ferraty and Philippe Vieu. *Nonparametric functional data analysis: theory and practice*. Springer Science & Business Media, 2006.

[6] James O Ramsay and Bernard W Silverman. *Applied functional data analysis: methods and case studies*. Springer, 2002.

[7] Mia Hubert, Peter J Rousseeuw, and Pieter Segaert. Multivariate functional outlier

detection. *Statistical Methods & Applications*,
24(2):177–202, 2015.

[8] Antonio Cuevas, Manuel Febrero, and Ricardo Fraiman. Robust estimation and classification for functional data via projection-based depth notions. *Computational Statistics*, 22(3):481–496, 2007.

[9] Guillaume Staerman, Pavlo Mozharovskyi, and Stéphan Clémençon. The area of the convex hull of sampled curves: a robust functional statistical depth measure. In *Proceedings of the 23rd International Conference on Artificial Intelligence and Statistics (AIS-TATS 2020)*, volume 108, pages 570–579, 2020.

[10] John W. Tukey. Mathematics and the pic-
turing of data. In *Proceedings of the Interna-
tional Congress of Mathematicians, Vancouver,
1975*, volume 2, pages 523–531, 1975.

[11] David L Donoho, Miriam Gasko, et al. Breakdown properties of location estimates based on halfspace depth and projected outlyingness. *The Annals of Statistics*, 20(4):1803–1827, 1992.

[12] Claudia Becker, Roland Fried, and Sonja Kuhnt. *Robustness and complex data structures: Festschrift in honour of Ursula Gather*. Springer Science & Business Media, 2014.

[13] Stanislav Nagy, Irène Gijbels, and Daniel Hlubinka. Depth-based recognition of shape outlying functions. *Journal of Computational and Graphical Statistics*, 26(4):883–893, 2017.

[14] Peter J Rousseeuw and Katrien Van Driessen. A fast algorithm for the minimum covariance determinant estimator. *Technometrics*, 41:212–223, 1999.

[15] Wolfgang Polonik. Minimum volume sets and generalized quantile processes. *Stochastic Processes and their Applications*, 69(1):1–24, 1997.

[16] Clayton Scott and Robert Nowak. Learning minimum volume sets. *Journal of Machine*
*Learning Research*, 7:665–704, 2006.

[17] Bernhard Schölkopf, John C. Platt, John Shawe-Taylor, Alex Smola, and Robert C. Williamson. Estimating the support of a high-dimensional distribution. *Neural Computation*, 13(7):1443–1471, 2001.

[18] Markus M. Breunig, Hans-Peter Kriegel, Raymond T. Ng, and Jörg Sander. LOF: Identifying density-based local outliers. In *Proceedings of the 2000 ACM SIGMOD International Conference on Management of Data*, volume 29, pages 93–104. ACM, 2000.

[19] F. T. Liu, K. M. Ting, and Z. Zhou. Isolation forest. In *Proceedings of the Eighth IEEE International Conference on Data Mining*, pages 413–422, 2008.

[20] Sahand Hariri, Matias Carrasco Kind, and Robert J. Brunner. Extended isolation forest. *IEEE Transactions on Knowledge and Data Engineering*, pages 1–1, 2019.

[21] Yijun Zuo and Robert Serfling. General notions of statistical depth function. *Ann. Statist.*, 28(2):461–482, 04 2000.

[22] Karl Mosler. Depth statistics. In Claudia Becker, Roland Fried, and Sonja Kuhnt, editors, *Robustness and Complex Data Structures: Festschrift in Honour of Ursula Gather*, pages 17–34. Springer. Berlin, 2013.

[23] James Kuelbs and Joel Zinn. Half-region depth for stochastic processes. *Journal of Multivariate Analysis*, 142:86 – 105, 2015.

[24] Alicia Nieto-Reyes and Heather Battey. A topologically valid definition of depth for functional data. *Statistical Science*, 2016.

[25] Irène Gijbels, Stanislav Nagy, et al. On a general definition of depth for functional data. *Statistical Science*, 32(4):630–639, 2017.

[26] Karl Mosler and Yulia Polyakova. General notions of depth for functional data. *arXiv:1208.1981*, 2018.

[27] Gerda Claeskens, Mia Hubert, Leen Slaets, and Kaveh Vakili. Multivariate functional

halfspace depth. *Journal of the American Statistical Association*, 109(505):411–423, 2014.

[28] Ricardo Fraiman and Graciela Muniz. Trimmed means for functional data. *Test*, 10(2):419–440, 2001.

[29] G. Brys, M. Hubert, and A. Struyf. A robust measure of skewness. *Journal of Computational and Graphical Statistics*, 13(4):996–1017, 2004.

[30] Tom Fawcett. An Introduction to ROC Analysis. *Letters in Pattern Recognition*, 27(8):861–874, 2006.

[31] Stephan Clémençon and Nicolas Vayatis. Nonparametric estimation of the precision-recall curve. In *ICML '09: Proceedings of the 26th Annual International Conference on Machine Learning*, page 185–192, 2009.

[32] Peter Segaert, Mia Hubert, Peter Rousseeuw, and J. Raymaekers. mrfdepth: Depth measures in multivariate, regression and functional settings. *R package version 1.0.11*, 2019.

[33] Nicholas Tarabelloni, Ana Arribas-Gil, Francesca Ieva, Anna Maria Paganoni, and Juan Romo. Roahd: Robust analysis of high dimensional data. *R package version 1.4.1*, 2018.

[34] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python. *Journal of Machine Learning Research*, 12:2825–2830, 2011.

[35] Rob J. Hyndman and Han Lin Shang. Rainbow plots, bagplots, and boxplots for functional data. *Journal of Computational and Graphical Statistics*, 19(1):29–45, 2010.

[36] Ying Sun and Marc G. Genton. Functional boxplots. *Journal of Computational and Graphical Statistics*, 20(2):316–334, 2011.
[37] Weiyi Xie, Sebastian Kurtek, Karthik Bharath, and Ying Sun. A geometric approach to visualization of variability in functional data. *Journal of the American Statistical Association*, 112(519):979–993, 2017.

[38] Ana Arribas-Gil and Juan Romo. Shape outlier detection and visualization for functional data: the outliergram. *Biostatistics*, 15(4):603–619, 2014.

[39] Peter J. Rousseeuw, Jakob Raymaekers, and Mia Hubert. A measure of directional outlyingness with applications to image data and video. *Journal of Computational and Graphical Statistics*, 27(2):345–359, 2018.

[40] Wenlin Dai and Marc Genton. Multivariate functional data visualization and outlier detection. *Journal of Computational and Graphical Statistics*, 03 2017.
# Appendix A Benchmarked datasets

In this part, we display in Figure B.1 the aeronautics and the rocks datasets.

# Appendix B Additional experiments on simulated anomalies

In this part, complementary experiments to the Section 3 are displayed. They are conducted with the same methodology but varying proportion of anomalies: 1% in Table B1, 2% in Table B2, 3% in Table B3 and 4% in Table B4.

<figure><img src="image_13.png" /><figcaption>Fig. B.1 The aeronautics and the rock datasets.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Rocks data</figcaption></figure>
<table><caption>**Table B1** Methods considered in performance comparison with the sensitivity ($p_c$) and the Area Under the Receiver Operating Characteristic (AUROC) for the four simulated models with 1% of added anomalies. Bold numbers correspond to the best result.</caption><thead><tr><th rowspan="2">Methods</th><th colspan="2">Isolated</th><th colspan="2">Magnitude I</th><th colspan="2">Magnitude II</th><th colspan="2">Shape</th></tr><tr><th>p<sub>c</sub></th><th>AUC</th><th>p<sub>c</sub></th><th>AUC</th><th>p<sub>c</sub></th><th>AUC</th><th>p<sub>c</sub></th><th>AUC</th></tr></thead><tbody><tr><td>FIF</td><td>0</td><td>0.20</td><td>1</td><td>1</td><td>0</td><td>0.32</td><td>0.06</td><td>0.98</td></tr><tr><td>fAO</td><td>0</td><td>0.44</td><td>1</td><td>1</td><td>0</td><td>0.54</td><td>0</td><td>0.67</td></tr><tr><td>fbd</td><td>0</td><td>0.44</td><td>1</td><td>1</td><td>0</td><td>0.54</td><td>0</td><td>0.68</td></tr><tr><td>fSDO</td><td>0</td><td>0.42</td><td>1</td><td>1</td><td>0</td><td>0.54</td><td>0</td><td>0.77</td></tr><tr><td>fT</td><td>0</td><td>0.43</td><td>1</td><td>1</td><td>0</td><td>0.44</td><td>0</td><td>0.72</td></tr><tr><td>ACH</td><td>0</td><td>0.62</td><td>1</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0.61</td></tr><tr><td>Outliergram</td><td>0</td><td>0.55</td><td>1</td><td>1</td><td>0</td><td>0.54</td><td>0</td><td>0.51</td></tr><tr><td>MS + IF</td><td>0</td><td>0.05</td><td>1</td><td>1</td><td>0</td><td>0.77</td><td>0</td><td>0.77</td></tr><tr><td>FOM (fSDO) + IF</td><td>0</td><td>0.27</td><td>1</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0.97</td></tr><tr><td>FOM (fAO) + IF</td><td>0</td><td>0.31</td><td>1</td><td>1</td><td>0.88</td><td>1</td><td>0</td><td>0.96</td></tr><tr><td>FPCA + IF</td><td>0</td><td>0.08</td><td>0</td><td>0.96</td><td>0</td><td>0.77</td><td>0</td><td>0.86</td></tr><tr><td>FPCA + LOF</td><td>0</td><td>0.31</td><td>0.18</td><td>0.86</td><td>0</td><td>0.42</td><td>0</td><td>0.67</td></tr><tr><td>FPCA + OC</td><td>0</td><td>0.03</td><td>0</td><td>0.93</td><td>0</td><td>0.78</td><td>0</td><td>0.88</td></tr></tbody></table>

<table><caption>**Table B2** Methods considered in performance comparison with the sensitivity ($p_c$) and the Area Under the Receiver Operating Characteristic (AUROC) for the four simulated models with 2% of added anomalies. Bold numbers correspond to the best result.</caption><thead><tr><th rowspan="2">Methods</th><th colspan="2">Isolated</th><th colspan="2">Magnitude I</th><th colspan="2">Magnitude II</th><th colspan="2">Shape</th></tr><tr><th>p<sub>c</sub></th><th>AUC</th><th>p<sub>c</sub></th><th>AUC</th><th>p<sub>c</sub></th><th>AUC</th><th>p<sub>c</sub></th><th>AUC</th></tr></thead><tbody><tr><td>FIF</td><td>0</td><td>0.23</td><td>0.97</td><td>1</td><td>0</td><td>0.32</td><td>0.26</td><td>0.98</td></tr><tr><td>fAO</td><td>0</td><td>0.44</td><td>1</td><td>1</td><td>0</td><td>0.54</td><td>0</td><td>0.67</td></tr><tr><td>fbd</td><td>0</td><td>0.44</td><td>1</td><td>1</td><td>0</td><td>0.54</td><td>0</td><td>0.68</td></tr><tr><td>fSDO</td><td>0</td><td>0.42</td><td>1</td><td>1</td><td>0</td><td>0.43</td><td>0</td><td>0.77</td></tr><tr><td>fT</td><td>0</td><td>0.43</td><td>1</td><td>1</td><td>0</td><td>0.44</td><td>0</td><td>0.72</td></tr><tr><td>ACH</td><td>0</td><td>0.62</td><td>1</td><td>0.83</td><td>0.94</td><td>1</td><td>0</td><td>0.61</td></tr><tr><td>Outliergram</td><td>0</td><td>0.55</td><td>1</td><td>1</td><td>0</td><td>0.54</td><td>0</td><td>0.49</td></tr><tr><td>MS + IF</td><td>0</td><td>0.05</td><td>1</td><td>1</td><td>0</td><td>0.77</td><td>0</td><td>0.77</td></tr><tr><td>FOM (fSDO) + IF</td><td>0</td><td>0.22</td><td>0.94</td><td>1</td><td>0.86</td><td>1</td><td>0</td><td>0.95</td></tr><tr><td>FOM (fAO) + IF</td><td>0</td><td>0.26</td><td>0.89</td><td>1</td><td>0.66</td><td>0.99</td><td>0</td><td>0.95</td></tr><tr><td>FPCA + IF</td><td>0</td><td>0.08</td><td>0</td><td>0.92</td><td>0</td><td>0.71</td><td>0</td><td>0.90</td></tr><tr><td>FPCA + LOF</td><td>0</td><td>0.35</td><td>0.09</td><td>0.81</td><td>0</td><td>0.48</td><td>0.09</td><td>0.77</td></tr><tr><td>FPCA + OC</td><td>0</td><td>0.03</td><td>0</td><td>0.94</td><td>0</td><td>0.79</td><td>0</td><td>0.91</td></tr></tbody></table>
<table><caption>**Table B3** Methods considered in performance comparison with the sensitivity ($p_c$) and the Area Under the Receiver Operating Characteristic (AUROC) for the four simulated models with 3% of added anomalies. Bold numbers correspond to the best result.</caption><thead><tr><th rowspan="2">Methods</th><th colspan="2">Isolated</th><th colspan="2">Magnitude I</th><th colspan="2">Magnitude II</th><th colspan="2">Shape</th></tr><tr><th>p<sub>c</sub></th><th>AUC</th><th>p<sub>c</sub></th><th>AUC</th><th>p<sub>c</sub></th><th>AUC</th><th>p<sub>c</sub></th><th>AUC</th></tr></thead><tbody><tr><td>FIF</td><td>0</td><td>0.21</td><td>0.98</td><td>1</td><td>0</td><td>0.33</td><td><strong>0.49</strong></td><td><strong>0.98</strong></td></tr><tr><td>fAO</td><td>0</td><td>0.44</td><td>1</td><td>1</td><td>0</td><td>0.54</td><td>0</td><td>0.67</td></tr><tr><td>fbd</td><td>0</td><td>0.44</td><td>1</td><td>1</td><td>0</td><td>0.54</td><td>0</td><td>0.68</td></tr><tr><td>fSDO</td><td>0</td><td>0.42</td><td>1</td><td>1</td><td>0</td><td>0.43</td><td>0</td><td>0.77</td></tr><tr><td>fT</td><td>0</td><td>0.43</td><td>1</td><td>1</td><td>0</td><td>0.44</td><td>0</td><td>0.72</td></tr><tr><td>ACH</td><td>0</td><td><strong>0.60</strong></td><td>1</td><td>0.85</td><td><strong>0.88</strong></td><td>1</td><td>0</td><td>0.60</td></tr><tr><td>Outliergram</td><td>0</td><td>0.55</td><td>1</td><td>1</td><td>0</td><td>0.54</td><td>0</td><td>0.49</td></tr><tr><td>MS + IF</td><td>0</td><td>0.05</td><td>1</td><td>1</td><td>0</td><td>0.75</td><td>0</td><td>0.77</td></tr><tr><td>FOM (fSDO) + IF</td><td>0</td><td>0.06</td><td>0.89</td><td>1</td><td>0.64</td><td>0.99</td><td>0</td><td>0.89</td></tr><tr><td>FOM (fAO) + IF</td><td>0</td><td>0.24</td><td>0.81</td><td>1</td><td>0.70</td><td>0.99</td><td>0</td><td>0.93</td></tr><tr><td>FPCA + IF</td><td>0</td><td>0.07</td><td>0</td><td>0.93</td><td>0</td><td>0.70</td><td>0</td><td>0.93</td></tr><tr><td>FPCA + LOF</td><td>0</td><td>0.39</td><td>0.16</td><td>0.87</td><td>0</td><td>0.57</td><td>0.17</td><td>0.70</td></tr><tr><td>FPCA + OC</td><td>0</td><td>0.04</td><td>0</td><td>0.94</td><td>0</td><td>0.78</td><td>0</td><td>0.93</td></tr></tbody></table>

<table><caption>**Table B4** Methods considered in performance comparison with the sensitivity ($p_c$) and the Area Under the Receiver Operating Characteristic (AUROC) for the four simulated models with 4% of added anomalies. Bold numbers correspond to the best result.</caption><thead><tr><th rowspan="2">Methods</th><th colspan="2">Isolated</th><th colspan="2">Magnitude I</th><th colspan="2">Magnitude II</th><th colspan="2">Shape</th></tr><tr><th>p<sub>c</sub></th><th>AUC</th><th>p<sub>c</sub></th><th>AUC</th><th>p<sub>c</sub></th><th>AUC</th><th>p<sub>c</sub></th><th>AUC</th></tr></thead><tbody><tr><td>FIF</td><td>0</td><td>0.23</td><td>0.94</td><td>1</td><td>0</td><td>0.34</td><td><strong>0.58</strong></td><td><strong>0.98</strong></td></tr><tr><td>fAO</td><td>0</td><td>0.44</td><td>1</td><td>1</td><td>0</td><td>0.54</td><td>0</td><td>0.67</td></tr><tr><td>fbd</td><td>0</td><td>0.44</td><td>1</td><td>1</td><td>0</td><td>0.54</td><td>0</td><td>0.68</td></tr><tr><td>fSDO</td><td>0</td><td>0.42</td><td>1</td><td>1</td><td>0</td><td>0.43</td><td>0</td><td>0.77</td></tr><tr><td>fT</td><td>0</td><td>0.43</td><td>1</td><td>1</td><td>0</td><td>0.44</td><td>0</td><td>0.72</td></tr><tr><td>ACH</td><td>0</td><td><strong>0.63</strong></td><td>1</td><td>0.85</td><td><strong>0.85</strong></td><td>1</td><td>0</td><td>0.60</td></tr><tr><td>Outliergram</td><td>0</td><td>0.55</td><td>1</td><td>1</td><td>0</td><td>0.54</td><td>0</td><td>0.46</td></tr><tr><td>MS + IF</td><td>0</td><td>0.05</td><td>1</td><td>1</td><td>0</td><td>0.74</td><td>0</td><td>0.77</td></tr><tr><td>FOM (fSDO) + IF</td><td>0</td><td>0</td><td>0.80</td><td>1</td><td>0.65</td><td>0.99</td><td>0</td><td>0.87</td></tr><tr><td>FOM (fAO) + IF</td><td>0</td><td>0.10</td><td>0.89</td><td>1</td><td>0.55</td><td>0.98</td><td>0</td><td>0.87</td></tr><tr><td>FPCA + IF</td><td>0</td><td>0.08</td><td>0</td><td>0.92</td><td>0</td><td>0.72</td><td>0.28</td><td>0.96</td></tr><tr><td>FPCA + LOF</td><td>0</td><td>0.43</td><td>0.26</td><td>0.82</td><td>0</td><td>0.59</td><td>0.18</td><td>0.70</td></tr><tr><td>FPCA + OC</td><td>0</td><td>0.03</td><td>0</td><td>0.93</td><td>0</td><td>0.78</td><td>0.08</td><td>0.95</td></tr></tbody></table>