Functional Anomaly Detection: a Benchmark Study

Introduction

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 contains ’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.

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.

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 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 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 ${\boldsymbol{X}}$ modelling the normal behavior of the system monitored takes its values in a Hilbert space, say the space $L^{2}()$ of square integrable functions on $$ 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 during the optimization, as pointed out in e.g. KuelbsZ15 , 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 nieto ; gijbels2017general 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 MoslerP18 , while integrating univariate depth functions over time is suggested in e.g. claeskens2014multivariate or hubert2015multivariate . This second approach is often preferred in practice, due to its simplicity and (often substantially) lower computational burden.

Integrated functional depths. Let $\text{D}^{(1)}(\cdot\mid\cdot)$ be a univariate depth measure, refer to e.g. ZuoS00 and Mosler13 for an account of depth statistics. A typical example is the Tukey depth, see tukey1975mathematics and also fraiman2001trimmed or claeskens2014multivariate . For an i.i.d. sample $\mathcal{D}_{n}=\{X_{1},\;\ldots,\;X_{n}\}$ of real valued observations, it is defined as:

where $\mathcal{C}_{n}(t)=\{{\boldsymbol{X}}_{1}(t),...,{\boldsymbol{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 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 ZuoS00 , which is given by:

where $\text{med}(X)$ and $\text{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 hubert2015multivariate is given by

where the asymmetric outlyingness AO is defined as:

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 BrysHS04 defined as follows:

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.

The Area of the Convex Hull (ACH) depth. The area of the convex hull (ACH) depth, recently introduced by staerman2019area , 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\leq J\leq n$, the ACH depth, denoted by $\text{D}_{ACH}({\boldsymbol{X}}\,\mid\,\mathcal{C}_{n})$ here, is defined as follows:

Functional Isolation Forest. Recently, staerman2019functional extended the popular Isolation Forest approach, see LiuTZ08 , to the functional framework. A Functional Isolation Forest (FIF) is formed by a collection of functional isolation trees (F-$i$trees), constructed each (by a sequence of random splits) from a subsample (of size $\psi$) of $\mathcal{C}_{n}$. The abnormality score of the observation ${\boldsymbol{X}}$ is then computed as a monotone decreasing transformation of the average depth of ${\boldsymbol{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 $\mathcal{C}_{n}$ with higher probability (and thus on less deep levels of the F-$i$trees) than those similar to the majority of the observed curves. The construction of the F-$i$trees 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 $\mathcal{C}_{n}$. Before each random univariate split, all observations of the F-$i$tree 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$:

where ${\boldsymbol{d}}\in\mathbf{D}$, ${\boldsymbol{X}}^{\prime}$ and $\|{\boldsymbol{X}}\|$ are respectively the first derivative and the $L_{2}$-norm of ${\boldsymbol{X}}$ in $L_{2}()$.

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.

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({\boldsymbol{x}})$. However, in unsupervised anomaly detection, the scoring function $s({\boldsymbol{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.

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., hubert2015multivariate . 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 ${\boldsymbol{X}}=({\boldsymbol{X}}(t))_{t\in}$ 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 ${\boldsymbol{Y}}$ (independently from the original dataset) that are next added to the selected above $1794$ ’normal’ observations ${\boldsymbol{X}}$. By $\mathcal{U}([a,b])$ is meant the uniform distribution on the interval $[a,b]$, while $\delta_{u}$ denotes the Dirac mass at point $u$.

${\boldsymbol{Y}}(t)\equiv u_{2}$, with $u_{2}\sim\mathcal{U}()$.

${\boldsymbol{Y}}(t)=\sin(2\pi u_{4}t)$, where $u_{4}\sim\mathcal{U}([0.2,2])$.

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

3 Naive Approaches - Sampled Curves Viewed as Multivariate Data

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

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 segaert ; the outliergram is available in the roahd R-package tarabelloni ; IF, LOF and OC are available in sklearn python library scikit-learn ; Magnitude-Shape plot (MS) is available in the scikit-fda python library https://github.com/GAA-UAM/scikit-fda; Functional Isolation Forest (FIF) and ACH open-source python codes are available under the following link https://github.com/GuillaumeStaermanML; 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).

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.

We analyze the complete Airbus dataset (i.e. including normal and abnormal accelerometer data) as well as data arising from the spectrophotometry 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.

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 hyndman10 , functional boxplots genton11 or amplitude and phase boxplot displays geomvisua , as well as those especially designed for anomaly detection such as the Outliergram 2014shape , the Functional Outlier Map (FOM) hubert2015multivariate ; FOM or the Magnitude-Shape plot (MS) MS . 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 $\mathcal{C}_{n}=\{{\boldsymbol{X}}_{1},...,{\boldsymbol{X}}_{n}\}$ observed on $\{t_{1},\ldots,t_{p}\}$, the MS plot is then the scatter plot of the points $(\text{MO}_{i},\text{VO}_{i})_{1\leq i\leq n}$ with the Mean Outlyingness (MO)

where the outyingness itself is usually a monotone decreasing transform of data depth: $\text{O}({\boldsymbol{X}}_{i}(t_{j}),\mathcal{C}_{n}(t_{j}))=\dfrac{1}{\text{D}({\boldsymbol{X}}_{i}(t_{j})\,\mid\,\mathcal{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 FOM . It is defined as a plot of points

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 auronautics data, and explain variations in the performance of different methods used in Section 4.2.

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.

Conclusion and Perspectives

Due to its unsupervised nature, the anomaly detection problem is extremely challenging. Because the learning algorithms used to build a 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 financial support from the initiative “Forschungspartnerschaften Mineralrohstoffe - ein strategischer Forschungsschwerpunkt der Geologischen Bundesanstalt”. The spectroscopic data of sedimentary material was provided by the Geological Survey of Austria.

References

Appendix A Benchmarked datasets

In this part, we display in Figure B.8 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 5, 2% in Table 6, 3% in Table 7 and 4% in Table 8.