# Quantum separability criteria based on realignment moments

Tinggui Zhang1,2† Naihuan Jing2,3 and Shao-Ming Fei1,4,5 ${1}$ School of
Mathematics and Statistics, Hainan Normal University, Haikou, 571158, China
$2$ College of Sciences, Shanghai University, Shanghai, 200444, China
$3$ Department of Mathematics, North Carolina State University
$4$ School of Mathematical Sciences, Capital Normal University, Beijing
100048, China
$5$ Max-Planck-Institute for Mathematics in the Sciences, Leipzig 04103,
Germany
† Correspondence to<EMAIL_ADDRESS>
###### Abstract

Abstract: Quantum entanglement is the core resource in quantum information
processing and quantum computing. It is an significant challenge to
effectively characterize the entanglement of quantum states. Recently, elegant
separability criterion is presented in [Phys. Rev. Lett. 125, 200501 (2020)]
by Elben _et al._ based on the first three partially transposed (PT) moments
of density matrices. Then in [Phys. Rev. Lett. 127, 060504 (2021)] Yu _et al._
proposed two general powerful criteria based on the PT moments. In this paper,
based on the realignment operations of matrices we propose entanglement
detection criteria in terms of such realignment moments. We show by detailed
example that the realignment moments can also be used to identify quantum
entanglement.

Keywords: Quantum entanglement; Separability criteria; Partially transposed
(PT) moments; Realignment moments

###### pacs:

03.67.-a, 02.20.Hj, 03.65.-w

##  I Introduction

Quantum entanglement is one of the most profound features of quantum mechanics
rpkh . It can be used to realize quantum tasks such as quantum communication
chgb ; akek ; cgcr , quantum simulation slio , quantum computation and so on
mail . How to characterize and quantify entanglement is one of the fundamental
problems in quantum information rpkh  and many-body quantum physics jeis .
There has been a lot of efforts to explore the quantum entanglement both
theoretically and experimentally rpkh ; oggt ; nfgv .

With the development of quantum information and quantum computing, the
research and development of medium-scale quantum devices containing dozens or
even more qubits is a hot topic jpre ; hszh ; mgon . For such devices, the
standard methods of tomography kvhr ; mpjr  and witnesses oggt ; nfgv  are no
longer feasible for gauging the performance in actual experiments mpjr ; yhua
. Such difficulties motivate people to construct direct entanglement
estimation scheme without any prior knowledge of the quantum states and
tomography. Recently, methods using locally randomized measurements to probe
entanglement have been put forward abcf ; tapb ; akno ; lkjw ; arhr ; snao ;
ypzl ; jlas ; hrjp ; xsog . Among them, many efforts have been devoted to
verify the positive partial transpose (PPT) condition aper  based on the
moments of the randomized measurements arhr ; ypzl ; jlas ; hrjp ; xsog .

Consider bipartite quantum states $\rho_{AB}$ in $H_{A}\otimes H_{B}$, where
$H_{A}$ and $H_{B}$ are $n$-dimensional Hilbert spaces. The well known PPT
criterion says that the partially transposed density matrix
$\rho_{AB}^{T_{A}}$ with respect to the subspace $H_{A}$ is positive
semidefinite if $\rho_{AB}$ is separable. Otherwise, $\rho_{AB}$ must be
entangled. Recently, it is found that the PPT condition can also be
illustrated by considering the so-called partial transpose moments (PT
moments),

$p_{k}\equiv Tr[(\rho_{AB}^{T_{A}})^{k}],~{}~{}~{}k=1,2,3\cdots,d,$

where $d=n^{2}$ is the dimension of the global system $H_{A}\otimes H_{B}$.
The PT moments can be experimentally measured through global (bilocal) random
unitary matrices jlas ; ypzl  or local randomized measurements arhr  based on
quantum shadow estimation hrjp . Obviously, $p_{1}=1$ for PPT states such that
$\rho_{AB}^{T_{A}}\geq 0$. PT moments been first studied in quantum field
theory in quantifying the correlations in many-body systems pjet . In arhr ,
the authors first show that the first three PT moments can be used to define a
simple yet powerful test for bipartite entanglement: $p_{3}\geq p_{2}^{2}$ if
$\rho_{AB}$ is a PPT state. More recently, the authors in xsog  obtain a
family of separability criteria: a necessary condition for $\rho_{AB}$ being a
separable state is that $B_{[(d-1)/2]}(\textbf{p})\geq 0$, where
$[B_{k}(\textbf{p})]_{ij}=p_{i+j+1}$ for $i,j=0,1,2,\cdots,k$,
$\textbf{p}=(p_{0},p_{1},p_{2},\cdots,p_{d})$ with $p_{0}=d$.

In this paper, instead of the PT moments, we introduce rearrangement moments
based on the rearrangement operations kclw ; orud . Similar to that given by
the PT moments, we present quantum separability criteria given by the
rearrangement moments. By detailed example, we show that the rearrangement
moments based criteria may detect entanglement.

##  II Separability criterion based on realignment moments

In the following we present separability criteria based on the realignment of
matrices, by using results from the matrix analysis racr ; racr1  and the
Hoelder’s inequality.

For a given $m\times n$ matrix $C=[c_{ij}]$ with entries $c_{ij}$, the vector
$vec(C)$ is defined by

$vec(C)=(c_{11},\cdots,c_{m1},c_{12},\cdots,c_{m2},\cdots,c_{1n},\cdots,c_{mn}).$

Let $Z$ be an $m\times m$ block matrix with $n\times n$ block matrices. The
realigned matrix $R(Z)$ of $Z$ is defined by

$\displaystyle R(Z)\equiv\left(\begin{array}[]{ccccccc}vec(Z_{1,1})\\\
\vdots\\\ vec(Z_{m,1})\\\ \vdots\\\ vec(Z_{1,m})\\\ \vdots\\\
vec(Z_{m,m})\end{array}\right).$ (8)

Let $\rho_{AB}$ be a bipartite state in $H_{A}\otimes H_{B}$. We define the
following quantities to be realignment moments of $\rho_{AB}$,

$r_{k}\equiv
Tr[R(\rho_{AB})R(\rho_{AB})^{{\dagger}}]^{\frac{k}{2}},~{}~{}~{}k=1,2,3\cdots,d.$

For convenience we set $r_{0}=d$. We have the following conclusion on quantum
separability in terms of the realignment moments.

###### Theorem 1

If the quantum state $\rho_{AB}$ is separable, then $r_{2}^{2}\leq r_{3}$.

Proof Let $M$ be a Hermitian $d\times d$ matrix with eigenvalue decomposition
$M=\sum_{i=1}^{d}\lambda_{i}|x_{i}\rangle\langle x_{i}|$, where $\lambda_{i}$
and $|x_{i}\rangle$ are the eigenvalues and eigenvetors of $M$, respectively.
The Schatten-p norm ($p\geq 1$) $\|M\|_{p}$ of $M$ is given by

$\|M\|_{p}\equiv(\sum_{i=1}^{d}|\lambda_{i}|^{p})^{1/p}=Tr(|M|^{p})^{1/p},$

where
$|M|=\sqrt{MM^{{\dagger}}}=\sqrt{M^{2}}=\sum_{i=1}^{d}|\lambda_{i}||x_{i}\rangle\langle
x_{i}|$. Set $M=(R(\rho_{AB})R(\rho_{AB})^{{\dagger}})^{\frac{1}{2}}$. By the
definition of realignment moments, we obtain
$\|M\|^{4}_{2}=\|R(\rho_{AB})R(\rho_{AB})^{{\dagger}}\|^{4}_{2}=r_{2}^{2}$,
$\|M\|_{1}=r_{1}$ and $\|M\|_{3}^{3}=r_{3}$. It has been shown in arhr  that
$\|M\|_{p}$ satisfies the following matrix norm inequality,
$\|M\|^{4}_{2}\leq\|M\|_{1}\|M\|_{3}^{3}$. That is for any state $\rho_{AB}$
we have $r_{2}^{2}\leq r_{1}r_{3}$.

Moreover, in kclw  it has been shown that if a bipartite state $\rho_{AB}$ is
separable, then the trace norm $||R(\rho_{AB})||$ of the realigned matrix
$R(\rho_{AB})$, namely, the sum of all the singular values of $R(\rho_{AB})$
satisfies $||R(\rho_{AB})||\equiv
Tr[R(\rho_{AB})R(\rho_{AB})^{{\dagger}}]^{1/2}\leq 1$. Hence, if $\rho_{AB}$
to be a separable state we have $r_{1}\leq 1$, we can draw the conclusion of
the theorem.

Denote $\textbf{r}=(r_{0},r_{1},r_{2},\cdots,r_{d})$. We can construct the
Hankel matrices $[H_{k}(\textbf{r})]_{ij}=r_{i+j}$ and
$[B_{k}(\textbf{r})]_{ij}=r_{i+j+1}$ for $i,j=0,1,2,\cdots,k$ xsog . A
stronger separability criterion involving higher order $r_{k}$ can be derived.

###### Theorem 2

If $\rho_{AB}$ is separable, then $\widehat{H_{k}}(\textbf{r})\geq 0$ for
$k=1,2\cdots,\lfloor\frac{d}{2}\rfloor$ and $\widehat{B_{l}}(\textbf{r})\geq
0$ for $l=1,2,\cdots,\lfloor\frac{d-1}{2}\rfloor$, where $\lfloor\cdot\rfloor$
is the integer function, $\widehat{H_{k}}$ and $\widehat{B_{l}}$ are matrices
obtained by replacing $r_{1}$ with $1$ in $[H_{k}(\textbf{r})]$ and
$[B_{k}(\textbf{r})]$, respectively.

Proof For any state $\rho_{AB}$ let
$Y=(R(\rho_{AB})R(\rho_{AB})^{{\dagger}})^{1/2}$. Denote
$\textbf{v}\equiv(v_{0},v_{1},\cdots,v_{\lfloor\frac{d}{2}\rfloor})=(I=Y^{0},Y,Y^{2},\cdots,Y^{\lfloor\frac{d}{2}\rfloor})$
and
$\textbf{u}\equiv(u_{0},u_{1},\cdots,u_{\lfloor\frac{d-1}{2}\rfloor})=(Y^{\frac{1}{2}},Y^{\frac{3}{2}},\cdots,Y^{\lfloor\frac{d-1}{2}\rfloor})$.
Under the Hilbert-Schmidt inner product of matrices, the elements of the Gram
matrices for v and u are given by

$\langle v_{i},v_{j}\rangle=Tr(Y^{i+j})=r_{i+j}$

and

$\langle
u_{i},u_{j}\rangle=Tr(Y^{i+\frac{1}{2}}Y^{j+\frac{1}{2}})=Tr(Y^{i+j+1})=r_{i+j+1},$

which are just the elements of the Hankel matrices $H_{k}(\textbf{r})$ and
$B_{k}(\textbf{r})$. As the Gram matrices are always positive semidefinite
racr1 , we have $H_{k}(\textbf{r})\geq 0$ for
$k=1,2\cdots,\lfloor\frac{d}{2}\rfloor$ and $B_{l}(\textbf{r})\geq 0$ for
$l=1,2,\cdots,\lfloor\frac{d-1}{2}\rfloor$ for any state. Now, if $\rho_{AB}$
to be a separable state, we have $r_{1}\leq 1$, then let all $r_{1}$s in
$\lfloor H_{k}(\textbf{r})\rfloor$ and $\lfloor B_{k}(\textbf{r})\rfloor$ are
replaced by $1$, we can draw the conclusion of the theorem.

In particular, for $l=1$ in Theorem 2 one gets
$B_{1}=\left(\begin{array}[]{cc}r_{1}&r_{2}\\\ r_{2}&r_{3}\\\
\end{array}\right)$, then $r_{1}$ replaced by $1$, we get
$\widehat{B_{1}}=\left(\begin{array}[]{cc}1&r_{2}\\\ r_{2}&r_{3}\\\
\end{array}\right)\geq 0$, which gives rise to the conclusion of Theorem 1.
Theorem 2 can be used as effective separability criteria. It implies that if
an inequalities $\widehat{H_{k}(\textbf{r})}\geq 0$ for
$k=1,2\cdots,\lfloor\frac{d}{2}\rfloor$ or $\widehat{B_{l}(\textbf{r})}\geq 0$
for $l=1,2,\cdots,\lfloor\frac{d-1}{2}\rfloor$ is violated, $\rho_{AB}$ must
be an entangled state.

In fact, based on matrix realignment, we can also define the following
moments,

$q_{k}\equiv
Tr[R(\rho_{AB}-\rho_{A}\otimes\rho_{B})R(\rho_{AB}-\rho_{A}\otimes\rho_{B})^{{\dagger}}]^{k/2}$

for $k=1,2,3\cdots,d.$ Set $\textbf{q}=(q_{0},q_{1},\cdots,q_{d})$, where
$q_{0}=d$. In a similar way we have the following results.

###### Theorem 3

For any state $\rho_{AB}$, we have $H_{k}(\textbf{q})\geq 0$ for
$k=1,2\cdots,\lfloor\frac{d}{2}\rfloor$ and $B_{l}(\textbf{q})\geq 0$ for
$l=1,2,\cdots,\lfloor\frac{d-1}{2}\rfloor$.

In zhcj ; ogog  the authors show that if a bipartite state $\rho_{AB}$ is
separable, then

$\|R(\rho_{AB}-\rho_{A}\otimes\rho_{B})\|\leq\sqrt{(1-Tr\rho_{A}^{2})(1-Tr\rho_{B}^{2})},$

where $\rho_{A}$ and $\rho_{B}$ are reduced matrices with respect to the
subsystems $A$ and $B$. Taking $l=1$ in Theorem 3 and using the above result,
we have

###### Corollary 1

If $\rho_{AB}$ is a separable state, then
$q_{2}^{2}\leq(1-Tr\rho_{A}^{2})(1-Tr\rho_{B}^{2})q_{3}$.

In the following we give a concrete example to illustrate our entanglement
criterion based on realignment moments. Let us consider the Werner state,

$\rho_{p}=p|\psi\rangle\langle\psi|+\frac{1-p}{4}I_{4},$

where $0\leq p\leq 1$ and
$|\psi\rangle=\frac{1}{\sqrt{2}}(|00\rangle+|11\rangle)$, $I_{4}$ is the
$4\times 4$ identity matrix. By straightforward calculation we get $f\equiv
r_{3}-r_{2}^{2}=\frac{1}{16}(6p^{3}+1-9p^{4}-6p^{2})$ as a function of $p$.
The inequality $r_{2}^{2}\leq r_{3}$ in Theorem 1 is violated for $p>0.44$,
see Fig. 1.

##  III Conclusions and discussions

We have introduced two kinds of moments based on the realignment operations.
Based on these realignment moments we have obtained entanglement criteria.
Although we have not yet found an example to show that the realignment moments
based criteria work better than other criteria like the PT moments based ones,
the approach may be used to to detect multipartite entanglement by using the
generalized realignment criteria lkcw ; bjtm . A scheme to experimentally
measure the realignment moments via measuring the purity $Tr(\rho_{AB}^{2})$
abcf ; tapb  or by global random unitary matrices jlas ; ypzl  or local
randomized measurements zlpz  would also be of significance.

Acknowledgments: This work is supported by the National Natural Science
Foundation of China under Grant Nos. 12126314, 12126351, 11861031, 12075159
and 12171044, Beijing Natural Science Foundation (Z190005), the Hainan
Provincial Natural Science Foundation of China under Grant No.121RC539, and
Academy for Multidisciplinary Studies, Capital Normal University. This project
is also supported by the specific research fund of The Innovation Platform for
Academicians of Hainan Province under Grant No. YSPTZX202215.

Data availability statement: Our manuscript has no associated data.

## References

  * (1) R. Horodecki, P. Horodecki, M. Horodecki and K. Horodecki, Quantum entanglement, Rev. Mod. Phys. 81, 865(2009).
  * (2) C.H. Bennett and G. Brassard, in Proceedings of IEEE International Conference on Computers, Systems, and Signal Processing, Bangalore, India(IEEE,1984), p.175.
  * (3) A.K. Ekert, Quantum cryptography based on Bell’s theorem, Phys. Rev. Lett. 67,661(1991).
  * (4) C.H. Bennett, G. Brassard,C.Crépeau, R. Jozsa, A. Peres and W.K. Wootters, Teleporting an unknown quantum state via dual classical and Einstein-Podolsky-Rosen channels, Phys. Rev. Lett. 70, 1895(1993).
  * (5) S. Lioyd, Universal Quantum simulators, Science, 273, 1073(1996).
  * (6) M.A. Nielsen and I.L. Chuang, Quantum Computation and Quantum Information: 10th Anniversary Edition, 10th ed.(Cambridge University Press, New York, NY, USA, 2011).
  * (7) J. Eisert, M. Cramer and M.B. Plenio, Colloquium: Area laws for the entanglement entropy, Rev. Mod. Phys. 82, 277(2010).
  * (8) O. Gühne and G.Tóth, Entanglement detection, Phys. Rep. 474, 1(2009).
  * (9) N. Friis, G. Vitagliano, M. Malik and M. Huber, Entanglement certification from theory to experiment, Nat. Rev. Phys. 1, 72(2019).
  * (10) J. Preskill, Quantum computing in the NISQ era and beyond, Quantum 2, 79 (2018).
  * (11) H.S. Zhong, et al. Quantum computational advantage using photons, Science, 370(6523), 1460 (2020); H.S. Zhong, et al. Phase-programmable Gaussian Boson sampling using stimulated squeezed light, Phys. Rev. Lett. 127, 180502 (2021).
  * (12) M. Gong, et al. Quantum walks on a programmable two-dimensional 62-qubit superconducting processor, Science, 372(6545), 948 (2021); Y. Wu, et al. Strong quantum computational advantage using a superconducting quantum processor, Phys. Rev. Lett. 127, 180501 (2021).
  * (13) K. Vogel and H. Risken, Determination of quasiprobability distributions in terms of probability distributions for the rotated quadrature phase, Phys. Rev. A 40, 2847 (1989).
  * (14) M. Paris and J. Rehacek, Quantum state estimation, Lect. Notes Phys. 649, 7 (2004).
  * (15) Y. Huang, Computing quantum discord is NP-complete, New J. Phys. 16, 033027 (2014).
  * (16) A. Elben, B. Vermersch, C.F. Roos and P. Zoller, Statistical correlations between locally randomized measurements: A toolbox for probing entanglement in many-body quantum states, Phys. Rev. A 99, 052323 (2019).
  * (17) T. Brydges, A. Elben, P. Jurcevic, B. Vermersch, C. Maier, B. P. Lanyon, P. Zoller, R. Blatt and C. F. Roos, Probing renyi entanglement entropy via randomized measurements, Science 364, 260 (2019).
  * (18) A. Ketterer, N. Wyderka and O. Güne, Characterizing multipartite entanglement with moments of random correlations, Phys. Rev. Lett. 122, 120505 (2019).
  * (19) L. Knips, J. Dziewior, W. Klobus, W. Laskowski, T. Paterek, P.J. Shadbolt, H. Weinfurter and J.D.A. Meinecke, Multipartite entanglement analysis from random correlations, npj Quantum Inf. 6, 51 (2020).
  * (20) S. Imai, N. Wyderka, A. Ketterer and O. Güne, Bound entanglement from randomized measurements, Phys. Rev. Lett. 126, 150501 (2021).
  * (21) A. Elben, R. Kueng, H.Y.R. Huang, R. van Bijnen, C. Kokail, M. Dalmonte, P. Calabrese, B. Kraus, J. Preskill, P. Zoller and B. Vermersch, Mixed-state entanglement from local randomized measurements, Phys. Rev. Lett. 125, 200501 (2020).
  * (22) Y. Zhou, P. Zeng and Z. Liu, Single-copies estimation of entanglement negativity, Phys. Rev. Lett. 125, 200502 (2020).
  * (23) J. Gray, L. Banchi, A. Bayat and S. Bose, Machine-learning-assisted many-body entanglement measurement, Phys. Rev. Lett. 121, 150503 (2018).
  * (24) H.Y. Huang, R. Kueng and J. Preskill, Predicting many properties of a quantum system from very few measurements, Nat. Phys. 16, 1050 (2020).
  * (25) X.D. Yu , S. Imai and O. Gühne, Optimal entanglement certification from moments of the partial transpose, Phys. Rev. Lett. 127, 060504 (2021).
  * (26) A. Peres, Separability criterion for density matrices, Phys. Rev. Lett. 77, 1413 (1996).
  * (27) P. Calabrese, J. Cardy and E. Tonni, Entanglement negativity in quantum field theory, Phys. Rev. Lett. 109, 130502 (2012).
  * (28) K. Chen and L.A. Wu, A matrix realignment method for recognizing entanglement, Quantum Information and Computation, 3(3), 193 (2003).
  * (29) O. Rudolph, Further Results on the cross norm criterion for separability, Quantum Information Processing, 4(3), 219 (2005).
  * (30) R.A. Horn and C.R. Johnson, Topics in Matrix Analysis (Cambridge University Press, New York, 1991).
  * (31) R.A. Horn and C.R. Johnson, Matrix Analysis (Cambridge University Press, Cambridge, England, 2012).
  * (32) N.P. Pitsianis, Ph.D. thesis, The Kronecker product in approximation and fast transform generation, Cornell University, New York, 1997.
  * (33) C.J. Zhang, Y.S. Zhang, S. Zhang, G.C. Guo, Entanglement detection beyond the cross-norm or realignment criterion, Phys. Rev. A 77, 060301(R) (2008)
  * (34) O. Gittsovich, O. Gühne, P. Hyllus, J. Eisert, Unifying several separability conditions using the covariance matrix criterion, Phys. Rev. A 78, 052319 (2008)
  * (35) P. Horodecki, Separability criterion and inseparable mixed states with positive partial transposition, Phys. Lett. A 232, 333 (1997)
  * (36) K. Chen and L.A. Wu, The generalized partial transposition criterion for separability of multipartite quantum states, Phys. Lett. A 306, 14, (2002).
  * (37) B. Jungnitsch, T. Moroder and O. Gühne, Taming multi-particle entanglement, Phys. Rev. Lett. 106, 190502(2011).
  * (38) Z. Liu, P. Zeng, Y. Zhou and M. Gu, Characterizing correlation within multipartite quantum systems via local randomized measurements, arXiv:2108.08045 [quant-ph](2021).


2021 citnum

[1]Guillaume Staerman

[1]LTCI, Télécom Paris, Institut Polytechnique de Paris, France

2]Department of Operations and Information Systems, University of Graz,
Austria

3] Airbus AI Research, Toulouse, France

# Functional Anomaly Detection: a Benchmark Study

<EMAIL_ADDRESS>Eric Adjakossa  Pavlo Mozharovskyi  Vera
Hofer  Jayant Sen Gupta  Stéphan Clémençon  *  [  [

(Received: date / Accepted: date)

###### 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\geq 1$)
is well-documented in the litterature 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. rousseeuw2018anomaly  or
staerman2019functional  and the references therein. Functional data are
available in the form of functions, images and shapes or more general objects
wang2016functional , 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. ramsay2004functional , ferraty2006nonparametric ,
wang2016functional . 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. ramsay2007applied  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.
ferraty2006nonparametric . 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 hubert2015multivariate
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. cuevas2007robust
or staerman2019area . The depth of a point within a data cloud was first
introduced in tukey1975mathematics , in order to define a notion of
multivariate median in $\mathbb{R}^{d}$ with $d\geq 2$. Afterwards, so as to
measure the depth of an arbitrary point $x$, donoho1992breakdown  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 preorder 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. becker2014robustness , with
applications to functional anomaly detection in particular, refer to e.g.,
hubert2015multivariate  or nagy2017depth .

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.

##  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. RousseeuwD99 . However, in
most practical situations, the probability distribution of the normal
observations is unknown and nonparametric approaches like Minimum Volume Set
estimation (MV-set) Polonik97 ; ScottN06 , One-Class Support Vector Machines
(OCSVM) SchoelkopfPSTSW05 , Local Outlier Factor (LOF BreunigKNS00 ) or
Isolation Forest (IF) LiuTZ08 ; HaririCKB19  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 tukey1975mathematics  (refer also to ZuoS00  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\geq 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\mid
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\mid P)\geq\tau\\}$ with $\tau>0$) are considered
as abnormal, see e.g. Mosler13 .

Properly reconstructed function
---
Poorly reconstructed function
Figure 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).

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.

Method | F1-Score | AP | AUC | $p_{c}$
---|---|---|---|---
Aeronautics dataset
IF | 0.71 | 0.69 | 0.49 | 0.27
LOF | 0.68 | 0.73 | 0.5 | 0.21
OCSVM | 0.78 | 0.77 | 0.63 | 0.45
Rocks dataset
IF | 0.971 | 0.978 | 0.528 | 0
LOF | 0.973 | 0.965 | 0.495 | 0.05
OCSVM | 0.971 | 0.965 | 0.686 | 0
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.

###  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 ${\boldsymbol{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 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:

$\text{D}^{(1)}_{T}(x\,\mid\,\mathcal{D}_{n})=\min\\{\hat{F}_{n}(x),1-\hat{F}_{n}(x^{-})\\}\,,$
(1)

where the càd-làg function
$\hat{F}_{n}(x)=(1/n)\sum_{i=1}^{n}\mathbb{I}\\{X_{i}\leq x\\}$ is the
empirical cdf based on the dataset $\mathcal{D}_{n}$, denoting by
$\mathbb{I}\\{\mathcal{E}\\}$ the indicator function of any event
$\mathcal{E}$. Now, given a sample of curves
$\mathcal{C}_{n}=\\{{\boldsymbol{X}}_{1},...,{\boldsymbol{X}}_{n}\\}$ composed
of independent copies of the functional observation ${\boldsymbol{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])$:

$\text{D}({\boldsymbol{X}}\mid\mathcal{C}_{n})=\int_{0}^{1}\text{D}^{(1)}\bigl{(}{\boldsymbol{X}}(t)\,\mid\,\mathcal{C}_{n}(t)\bigr{)}dt,$
(2)

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 $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 ZuoS00 , which is given by:

$\text{D}^{(1)}_{P}(x\,\mid\,\mathcal{D}_{n})=\left(1+\frac{\lvert
x-\text{med}(X)\rvert}{\text{MAD}(X)}\right)^{-1}\,,$

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

$\text{D}^{(1)}_{AP}(x\,\mid\,\mathcal{D}_{n})=\frac{1}{1+\text{AO}(x\,\mid\,\mathcal{D}_{n})}\,,$

where the asymmetric outlyingness AO is defined as:

$\text{AO}(x\,\mid\,\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

$\displaystyle w_{1}(X)$
$\displaystyle=\hat{q}_{X}(0.25)-IQR(X)1.5e^{-4\text{MC}(X)}$ $\displaystyle
w_{2}(X)$ $\displaystyle=\hat{q}_{X}(0.75)+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 BrysHS04  defined as follows:

$\text{med}\left(\frac{\bigl{(}X_{j}+X_{i}-2~{}\text{med}(X)\bigr{)}}{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 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).

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:

${\small\frac{1}{\displaystyle\binom{n}{J}}\sum_{1\leq i_{1}\leq\dots\leq
i_{j}\leq
n}\frac{\lambda(\text{conv}(\text{graph}(\\{{\boldsymbol{X}}_{i_{1}},\dots,{\boldsymbol{X}}_{i_{j}}\\})))}{\lambda(\text{conv}(\text{graph}(\\{{\boldsymbol{X}}_{i_{1}},\dots,{\boldsymbol{X}}_{i_{j}}\\}\cup\\{{\boldsymbol{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}(\\{{\boldsymbol{X}}_{1},\dots,{\boldsymbol{X}}_{n}\\})=\bigcup_{i=1}^{n}\text{graph}({\boldsymbol{X}}_{i})$
with the notation
$\text{graph}({\boldsymbol{X}}_{i})=\\{(t,y):y={\boldsymbol{X}}_{i}(t),t\in[0,1]\\}$.
In staerman2019area  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 preorder induced by the resulting depth
function. One may refer to staerman2019area  for a detailed description of
this approach.

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[0,1]$:

$\langle{\boldsymbol{X}},{\boldsymbol{d}}\rangle:=\alpha\times\frac{\langle{\boldsymbol{X}},{\boldsymbol{d}}\rangle_{L_{2}}}{\|{\boldsymbol{X}}\|\|{\boldsymbol{d}}\|}+(1-\alpha)\times\frac{\langle{\boldsymbol{X}}^{\prime},{\boldsymbol{d}}^{\prime}\rangle_{L_{2}}}{\|{\boldsymbol{X}}^{\prime}\|\|{\boldsymbol{d}}^{\prime}\|},$

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}([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
${\boldsymbol{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}\rightarrow\mathbb{R}\cup\\{+\infty\\}$ in an unsupervised
manner (without observing $Y$), so that, ideally, the larger
$s({\boldsymbol{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({\boldsymbol{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. ${\boldsymbol{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 vs the
true positive rate, namely
$t\in\mathbb{R}\mapsto(\mathbb{P}\\{s({\boldsymbol{X}})\geq t\mid
Y=-1\\},\;\mathbb{P}\\{s({\boldsymbol{X}})\geq t\mid Y=+1\\})$, referred to as
the ROC curve (standing for Receiver Operator Characteristic curve), see e.g.
Fawcett06 . 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({\boldsymbol{X}}):=(dF_{+}/dF_{-})({\boldsymbol{X}})$, denoting by
$F_{\sigma}$ the conditional distribution of ${\boldsymbol{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:
$\mathrm{AUC}(s)=\mathbb{P}\\{s({\boldsymbol{X}})>s({\boldsymbol{X}}^{\prime})\mid
Y=+1,\;Y^{\prime}=-1\\}+\mathbb{P}\\{s({\boldsymbol{X}})=s({\boldsymbol{X}}^{\prime})\mid
Y=+1,\;Y^{\prime}=-1\\}/2$, where $({\boldsymbol{X}}^{\prime},Y^{\prime})$
denotes an independent copy of the pair $({\boldsymbol{X}},Y)$.

PR analysis. Alternatively, one may evalute the accuracy of any score by
plotting the precision-recall (shortly PR) curve, namely
$t\mapsto(\mathbb{P}\\{s({\boldsymbol{X}})\geq t\mid
Y=+1\\},\;\mathbb{P}\\{Y=+1\mid s({\boldsymbol{X}})\geq 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. CV11 . 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({\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.

###  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.,
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[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 ${\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$.

Model 1 (Isolated Anomalies)
${\boldsymbol{Y}}(t)=\varepsilon u_{1}\mathbb{I}\\{t=\tau\\}$, where
$u_{1}\sim\mathcal{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)

${\boldsymbol{Y}}(t)\equiv u_{2}$, with $u_{2}\sim\mathcal{U}([-12,-15])$.

Model 3 (Magnitude Anomalies II)

${\boldsymbol{Y}}(t)=u_{3}\mathbb{I}(t\in I)$, where
$u_{3}\sim\mathcal{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)

${\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.

Anomaly type |  | IF | LOF | OCSVM
---|---|---|---|---
Isolated | $p_{c}$ | 0 | 0 | 0
AUC | 0.41 | 0.25 | 0.44
Magnitude I | $p_{c}$ | 1 | 0.48 | 1
AUROC | 1 | 0.97 | 1
Magnitude II | $p_{c}$ | 0 | 0 | 0
AUC | 0.54 | 0.03 | 0.7
Shape | $p_{c}$ | 0 | 0.48 | 0
AUC | 0.67 | 0.97 | 0.67
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. Isolated |  Magnitude I
---|---
|
Magnitude II |  Shape
|
Figure 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. Isolated |  Magnitude I
---|---
|
Magnitude II |  Shape
|
Figure 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. | Isolated | Magnitude I | Magnitude II | Shape
---|---|---|---|---
Methods | $p_{c}$ | AUC | $p_{c}$ | AUC | $p_{c}$ | AUC | $p_{c}$ | AUC
FIF | 0 | 0.21 | 0.93 | 0.99 | 0 | 0.30 | 0.63 | 0.98
fAO | 0 | 0.44 | 1 | 1 | 0 | 0.54 | 0 | 0.66
fbd | 0 | 0.44 | 1 | 1 | 0 | 0.54 | 0 | 0.68
fSDO | 0 | 0.42 | 1 | 1 | 0 | 0.43 | 0 | 0.77
fT | 0 | 0.43 | 1 | 1 | 0 | 0.44 | 0 | 0.71
ACH | 0 | 0.63 | 0.48 | 0.97 | 0.80 | 0.99 | 0 | 0.56
Outliergram | 0 | 0.55 | 1 | 1 | 0 | 0.54 | 0 | 0.47
MS + IF | 0 | 0.05 | 0.66 | 0.98 | 0 | 0.70 | 0.33 | 0.74
FOM (fSDO) + IF | 0 | 0 | 0.85 | 0.99 | 0.64 | 0.99 | 0.06 | 0.80
FOM (fAO) + IF | 0 | 0.14 | 0.81 | 0.99 | 0.55 | 0.98 | 0.02 | 0.87
FPCA + IF | 0 | 0.11 | 0 | 0.91 | 0 | 0.71 | 0.45 | 0.97
FPCA + LOF | 0 | 0.5 | 0 | 0.16 | 0 | 0.38 | 0 | 0.79
FPCA + OC | 0 | 0.04 | 0 | 0.93 | 0 | 0.77 | 0.3 | 0.96
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.  Isolated | Magnitude I
---|---
|
Magnitude II | Shape
|
Figure 5: ROC curves of the four best methods for the simulation study.

###  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 claeskens2014multivariate .
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 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 111https://github.com/GAA-UAM/scikit-fda; Functional Isolation Forest
(FIF) and ACH open-source python codes are available under the following link
222https://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).

##  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.

FPCA |  MS-plot
---|---
|
FOM-fAO |  FOM-fSDO
|
Figure 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. | Airbus data | Rocks data
---|---|---
Method | F1-Score | AP | AUC | $p_{c}$ | F1-Score | AP | AUC | $p_{c}$
FIF | 0.81 | 0.88 | 0.76 | 0.52 | 0.974 | 0.992 | 0.772 | 0.10
fAO | 0.78 | 0.77 | 0.63 | 0.45 | 0.989 | 0.991 | 0.833 | 0.58
fbd | 0.78 | 0.77 | 0.63 | 0.45 | 0.977 | 0.988 | 0.751 | 0.19
fSDO | 0.78 | 0.77 | 0.63 | 0.45 | 0.972 | 0.980 | 0.555 | 0
fT | 0.78 | 0.77 | 0.63 | 0.45 | 0.984 | 0.989 | 0.780 | 0.43
ACH | 0.77 | 0.77 | 0.62 | 0.44 | 0.972 | 0.961 | 0.280 | 0
Outliergram | 0.71 | 0.76 | 0.55 | 0.28 | 0.974 | 0.981 | 0.66 | 0.03
MS + IF | 0.80 | 0.76 | 0.64 | 0.51 | 0.972 | 0.981 | 0.601 | 0
FOM (fSDO) + IF | 0.81 | 0.76 | 0.66 | 0.53 | 0.972 | 0.978 | 0.530 | 0.02
FOM (fAO) + IF | 0.80 | 0.77 | 0.65 | 0.51 | 0.984 | 0.991 | 0.804 | 0.448
FPCA + IF | 0.81 | 0.80 | 0.70 | 0.53 | 0.972 | 0.971 | 0.446 | 0
FPCA + LOF | 0.72 | 0.73 | 0.52 | 0.31 | 0.972 | 0.969 | 0.445 | 0.02
FPCA + OC | 0.81 | 0.79 | 0.70 | 0.54 | 0.972 | 0.971 | 0.463 | 0
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.

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.

###  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 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)

$\text{MO}_{i}=\frac{1}{p}\sum_{j=1}^{p}\;\text{O}({\boldsymbol{X}}_{i}(t_{j})\,\mid\,\mathcal{C}_{n}(t_{j})),$

and the Variational Outlyingness (VO)

$\text{VO}_{i}=\frac{1}{p}\sum_{j=1}^{p}\;\left(\text{O}({\boldsymbol{X}}(t_{j})\,\mid\,\mathcal{C}_{n}(t_{j}))-\text{MO}_{i}\right)^{2},$

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

$\Bigl{(}\text{MO}_{i},\frac{\sqrt{\text{VO}_{i}}}{1+\text{MO}_{i}}\Bigr{)}_{1\leq
i\leq 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
auronautics 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.

---
Figure 7: ROC curves of the four best methods for the aeronautics (top) and
rocks (bottom) data.

##  5 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

  * [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 Application, 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 23nd International Conference on Artificial Intelligence and Statistics (AISTATS 2020), volume 108, pages 570–579, 2020.
  * [10] John W. Tukey.  Mathematics and the picturing of data.  In Proceedings of the International 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] Wolgang 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.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.

Aeronautics data | Rocks data
---|---
|
Figure B.8: The aeronautics and the rock datasets. | Isolated | Magnitude I | Magnitude II |  Shape
---|---|---|---|---
Methods | $p_{c}$ | AUC | $p_{c}$ | AUC | $p_{c}$ | AUC | $p_{c}$ | AUC
FIF | 0 | 0.20 | 1 | 1 | 0 | 0.32 | 0.06 | 0.98
fAO | 0 | 0.44 | 1 | 1 | 0 | 0.54 | 0 | 0.67
fbd | 0 | 0.44 | 1 | 1 | 0 | 0.54 | 0 | 0.68
fSDO | 0 | 0.42 | 1 | 1 | 0 | 0.54 | 0 | 0.77
fT | 0 | 0.43 | 1 | 1 | 0 | 0.44 | 0 | 0.72
ACH | 0 | 0.62 | 1 | 1 | 1 | 1 | 0 | 0.61
Outliergram | 0 | 0.55 | 1 | 1 | 0 | 0.54 | 0 | 0.51
MS + IF | 0 | 0.05 | 1 | 1 | 0 | 0.77 | 0 | 0.77
FOM (fSDO) + IF | 0 | 0.27 | 1 | 1 | 1 | 1 | 0 | 0.97
FOM (fAO) + IF | 0 | 0.31 | 1 | 1 | 0.88 | 1 | 0 | 0.96
FPCA + IF | 0 | 0.08 | 0 | 0.96 | 0 | 0.77 | 0 | 0.86
FPCA + LOF | 0 | 0.31 | 0.18 | 0.86 | 0 | 0.42 | 0 | 0.67
FPCA + OC | 0 | 0.03 | 0 | 0.93 | 0 | 0.78 | 0 | 0.88
Table 5: 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.  | Isolated | Magnitude I | Magnitude II |  Shape
---|---|---|---|---
Methods | $p_{c}$ | AUC | $p_{c}$ | AUC | $p_{c}$ | AUC | $p_{c}$ | AUC
FIF | 0 | 0.23 | 0.97 | 1 | 0 | 0.32 | 0.26 | 0.98
fAO | 0 | 0.44 | 1 | 1 | 0 | 0.54 | 0 | 0.67
fbd | 0 | 0.44 | 1 | 1 | 0 | 0.54 | 0 | 0.68
fSDO | 0 | 0.42 | 1 | 1 | 0 | 0.43 | 0 | 0.77
fT | 0 | 0.43 | 1 | 1 | 0 | 0.44 | 0 | 0.72
ACH | 0 | 0.62 | 1 | 0.83 | 0.94 | 1 | 0 | 0.61
Outliergram | 0 | 0.55 | 1 | 1 | 0 | 0.54 | 0 | 0.49
MS + IF | 0 | 0.05 | 1 | 1 | 0 | 0.77 | 0 | 0.77
FOM (fSDO) + IF | 0 | 0.22 | 0.94 | 1 | 0.86 | 1 | 0 | 0.95
FOM (fAO) + IF | 0 | 0.26 | 0.89 | 1 | 0.66 | 0.99 | 0 | 0.95
FPCA + IF | 0 | 0.08 | 0 | 0.92 | 0 | 0.71 | 0 | 0.90
FPCA + LOF | 0 | 0.35 | 0.09 | 0.81 | 0 | 0.48 | 0.09 | 0.77
FPCA + OC | 0 | 0.03 | 0 | 0.94 | 0 | 0.79 | 0 | 0.91
Table 6: 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.  | Isolated | Magnitude I | Magnitude II |  Shape
---|---|---|---|---
Methods | $p_{c}$ | AUC | $p_{c}$ | AUC | $p_{c}$ | AUC | $p_{c}$ | AUC
FIF | 0 | 0.21 | 0.98 | 1 | 0 | 0.33 | 0.49 | 0.98
fAO | 0 | 0.44 | 1 | 1 | 0 | 0.54 | 0 | 0.67
fbd | 0 | 0.44 | 1 | 1 | 0 | 0.54 | 0 | 0.68
fSDO | 0 | 0.42 | 1 | 1 | 0 | 0.43 | 0 | 0.77
fT | 0 | 0.43 | 1 | 1 | 0 | 0.44 | 0 | 0.72
ACH | 0 | 0.60 | 1 | 0.85 | 0.88 | 1 | 0 | 0.60
Outliergram | 0 | 0.55 | 1 | 1 | 0 | 0.54 | 0 | 0.49
MS + IF | 0 | 0.05 | 1 | 1 | 0 | 0.75 | 0 | 0.77
FOM (fSDO) + IF | 0 | 0.06 | 0.89 | 1 | 0.64 | 0.99 | 0 | 0.89
FOM (fAO) + IF | 0 | 0.24 | 0.81 | 1 | 0.70 | 0.99 | 0 | 0.93
FPCA + IF | 0 | 0.07 | 0 | 0.93 | 0 | 0.70 | 0 | 0.93
FPCA + LOF | 0 | 0.39 | 0.16 | 0.87 | 0 | 0.57 | 0.17 | 0.70
FPCA + OC | 0 | 0.04 | 0 | 0.94 | 0 | 0.78 | 0 | 0.93
Table 7: 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.  | Isolated | Magnitude I | Magnitude II |  Shape
---|---|---|---|---
Methods | $p_{c}$ | AUC | $p_{c}$ | AUC | $p_{c}$ | AUC | $p_{c}$ | AUC
FIF | 0 | 0.23 | 0.94 | 1 | 0 | 0.34 | 0.58 | 0.98
fAO | 0 | 0.44 | 1 | 1 | 0 | 0.54 | 0 | 0.67
fbd | 0 | 0.44 | 1 | 1 | 0 | 0.54 | 0 | 0.68
fSDO | 0 | 0.42 | 1 | 1 | 0 | 0.43 | 0 | 0.77
fT | 0 | 0.43 | 1 | 1 | 0 | 0.44 | 0 | 0.72
ACH | 0 | 0.63 | 1 | 0.85 | 0.85 | 1 | 0 | 0.60
Outliergram | 0 | 0.55 | 1 | 1 | 0 | 0.54 | 0 | 0.46
MS + IF | 0 | 0.05 | 1 | 1 | 0 | 0.74 | 0 | 0.77
FOM (fSDO) + IF | 0 | 0 | 0.80 | 1 | 0.65 | 0.99 | 0 | 0.87
FOM (fAO) + IF | 0 | 0.10 | 0.89 | 1 | 0.55 | 0.98 | 0 | 0.87
FPCA + IF | 0 | 0.08 | 0 | 0.92 | 0 | 0.72 | 0.28 | 0.96
FPCA + LOF | 0 | 0.43 | 0.26 | 0.82 | 0 | 0.59 | 0.18 | 0.70
FPCA + OC | 0 | 0.03 | 0 | 0.93 | 0 | 0.78 | 0.08 | 0.95
Table 8: 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.