Elastic depths for detecting shape anomalies in functional data

TL;DR

This paper introduces elastic depths, a new family of depth measures that effectively detect shape anomalies in functional data by leveraging elastic distances in amplitude and phase spaces, outperforming existing methods.

Contribution

The paper presents elastic depths, a novel approach for shape anomaly detection in functional data, with theoretical advantages and practical implementation, including an R-package.

Findings

Elastic depths improve shape anomaly detection accuracy.

Compared to existing methods, elastic depths show superior performance.

Demonstrated effectiveness on hurricane trajectory data.

Abstract

We propose a new family of depth measures called the elastic depths that can be used to greatly improve shape anomaly detection in functional data. Shape anomalies are functions that have considerably different geometric forms or features from the rest of the data. Identifying them is generally more difficult than identifying magnitude anomalies because shape anomalies are often not distinguishable from the bulk of the data with visualization methods. The proposed elastic depths use the recently developed elastic distances to directly measure the centrality of functions in the amplitude and phase spaces. Measuring shape outlyingness in these spaces provides a rigorous quantification of shape, which gives the elastic depths a strong theoretical and practical advantage over other methods in detecting shape anomalies. A simple boxplot and thresholding method is introduced to identify shape anomalies using the elastic depths. We assess the elastic depth's detection skill on simulated shape outlier scenarios and compare them against popular shape anomaly detectors. Finally, we use hurricane trajectories to demonstrate the elastic depth methodology on manifold valued functional data. Supplementary materials, including additional simulations, data examples, and an R-package are available online.

Tables

Table 1. Table 1: Average rank of the single outlier by detection method and outlier model. Lower ranks (closer to 1) indicate the detection method ranked the outlier as more outlying. Bold font indicates the top results in each outlier model.

	Method	Model 1	Model 2	Model 3	Model 4	Model 5	Model 6	Model 7
1	ED-A	1.000	1.002	1.002	1.016	1.000	1.001	8.808
2	ED-P	62.832	1.027	1.449	2.012	2.831	11.151	2.689
3	TVD	1.000	26.341	4.453	6.228	1.000	1.226	2.355
4	DIR	1.024	40.764	4.678	33.504	4.380	1.011	10.082
5	MS	1.000	44.306	2.732	36.242	2.370	1.000	9.412
6	GEOM	1.000	99.995	93.343	30.776	1.002	7.590	92.303

Table 2. Table 2: Average rank of the single outlier by detection method and outlier model. Lower ranks (closer to 1) indicate the detection method ranked the outlier as more outlying. Bold font indicates the top results in each outlier model.

	Method	Model 1	Model 2	Model 3	Model 4	Model 5	Model 6	Model 7
1	ED-A	1.000	1.002	1.002	1.016	1.000	1.001	8.808
2	ST-T1	1.000	20.125	1.226	15.361	1.058	1.000	1.979
3	TVD	1.000	26.341	4.453	6.228	1.000	1.226	2.355
4	ST-T2	1.285	13.582	1.096	15.700	1.712	10.909	1.591
5	rFUNTA	1.430	23.192	9.278	7.564	1.000	2.044	1.802
6	FOM	1.976	20.389	3.577	16.464	2.035	7.014	3.014
7	ED-P	62.832	1.027	1.449	2.012	2.831	11.151	2.689
8	DIR	1.024	40.764	4.678	33.504	4.380	1.011	10.082
9	MS	1.000	44.306	2.732	36.242	2.370	1.000	9.412
10	OG	9.887	35.803	13.180	29.691	12.083	10.004	17.344
11	FDJ	9.576	40.953	15.899	38.746	14.959	11.190	12.002
12	IDJ	33.140	40.763	41.528	38.871	40.575	38.972	42.730
13	GEOM	1.000	99.995	93.343	30.776	1.002	7.590	92.303
14	ST-D1	51.031	50.828	51.286	51.779	52.746	52.656	53.573

Equations

$$q_f(t)=\frac{f^{\prime }(t)}{\sqrt{|f^{\prime }(t)|}}.$$

$$||q_f-q_g|{|}_2=\sqrt{{\int }_0^1|q_f(t)-q_g(t){|}^{2}dt},$$

$$||q_{f\circ \gamma }-q_{g\circ \gamma }|{|}_2=||q_f-q_g|{|}_2,$$

$${\gamma }^{\ast }=\arg \underset{\gamma \in \Gamma }{\inf }||q_f-q_{g\circ \gamma }|{|}_2,$$

$$d_a(f,g)=\underset{\gamma \in \Gamma }{\inf }||q_f-q_{g\circ \gamma }|{|}_2,$$

$$F_M=\{f:[0,1]\mapsto M,f\text{ is differentiable and M is a Riemannian manifold}\}.$$

$$O_a(f,P)$$

$$O_p(f,P)$$

$$D_a(f,P)$$

$$D_p(f,P)$$

$$D_a(f\circ \gamma ,P)=D_a(f,P),$$

$$O_{a,n}(f,P_n)$$

$$O_{p,n}(f,P_n)$$

$$D_{a,n}(f,P_n)$$

$$D_{p,n}(f,P_n)$$

$$\underset{n\to \infty }{\lim }\underset{f\in F_M}{\sup }|D_{a,n}(f,P_n)-D_a(f,P)|$$

$$\underset{n\to \infty }{\lim }\underset{f\in F_M}{\sup }|D_{p,n}(f,P_n)-D_p(f,P)|$$

$$F_1$$

$$\Gamma =\{\gamma :[0,1]\mapsto [0,1]|\gamma (0)=0\text{, }\gamma (1)=1\text{, }\gamma \text{ is diffeomorphic}\}.$$

$$f(t)\mapsto f(t)/||f^{\prime }||\forall t\in [0,1],$$

$$O(f)=\{Of(t):\forall t\in [0,1]\},$$

$$q_f(t)=\frac{f^{\prime }(t)}{\sqrt{|f^{\prime }(t)|}},$$

$$d(q_f,q_g)=\arccos {\int }_0^1⟨q_f(t),q_g(t)⟩dt,$$

$$d_a(f,g)=\underset{\gamma \in \Gamma ,O\in SO(n)}{\inf }\arccos {\int }_0^1⟨q_f(t),q_{O(g\circ \gamma }(t)⟩dt,$$

$$f^{\prime }(t{)}_{f(t)\mapsto c}=f^{\prime }(t)-2⟨f^{\prime }(t),c⟩\frac{f(t)+c}{||f(t)+c|{|}^2}.$$

$$h_f(t)=\frac{f^{\prime }(t{)}_{f(t)\mapsto c}}{\sqrt{|f^{\prime }(t)|}},$$

$$d_a(f,g)=\underset{\gamma \in \Gamma }{\inf }||h_f-h_{g\circ \gamma }|{|}_2,,$$

$$d_p(f_1,f_2)=\arccos {\int }_0^1⟨\sqrt{I^{\prime }(t)},\sqrt{{\gamma }^{\prime }(t)}⟩dt,$$

$$D_a(f\circ \gamma ,P)$$

$$=\text{median}\{d_a(f\circ \gamma ,X):X\in F_M\}$$

$$=\text{median}\{d_a(f,X):X\in F_M\}$$

$$=(1+O(f,P){)}^{-1}$$

$$=D_a(f,P)$$

$$s$$

$$=\arg \underset{f\in F_M}{\min }{O}_a(f,P)$$

Full text

Keywords: Anomaly detection; Data depth; Functional data; Shape analysis

1 Introduction

As data collection methods rapidly advance, functional data and functional data analysis (FDA) have become more prevalent. Functional data refers to data collected continuously across a compact domain, such as a fixed length of time or region of space, and where an observation is an entire curve or surface over the domain, rather than a single value. Examples of functional data include growth rate curves, electrocardiogram (ECG) data, temperature profiles, imaging data containing geometric shapes, and hurricane trajectories (See Figure 1).

As with traditional data analysis methods, it is critical to perform exploratory data analysis with functional data. Exploratory analysis can reveal significant trends or anomalies, which could bias post-processing analysis such as model fitting. Functional anomalies are of particular interest because of the adverse effects they can have on statistical models. Functional anomalies can also be interesting in their own right and can even be the primary focus of study.

The hurricane trajectories in Figure 1 present a challenge to current functional shape anomaly detectors. Due to the vast distances hurricanes travel, their trajectories resemble paths along the surface of a sphere. Existing methods do not handle spherical valued data, or more generally manifold valued data, so they have to approximate these paths with two-dimensional trajectories. This approximation distorts the distances between curves and consequently has a strong influence on the detection of anomalies. Furthermore, exploratory analysis shows that these trajectories exhibit substantial phase variability that could affect shape outlier detection (Srivastava et al.,, 2011).

In the functional data setting, identifying functional anomalies, i.e. identifying an entire function as an outlier, is not as straightforward as identifying univariate outliers with visualization methods. By definition, a functional anomaly is a function that is significantly more “extreme” in its characteristics than the rest of the functional data. Generally, functional anomalies are categorized into two types: magnitude and shape anomalies (Dai and Genton,, 2018). Magnitude anomalies are functions that clearly lie outside of the range of all other functions and are usually detected through data visualization methods (e.g., Hyndman and Shang,, 2010; Sun and Genton,, 2011; Myllymäki et al.,, 2017). On the other hand, shape anomalies take on a different shape or pattern than the rest of the data. They are more challenging to identify with visualization methods because they can lie hidden amongst the rest of the functions (Arribas-Gil and Romo,, 2014). Examples of shape anomalies include trajectories with more or less curvature, trajectories with more or fewer oscillations, or trajectories sampled from a process having a different mean function than the rest of the data. These are only a few examples; the possible ways a trajectory can be shape outlying are innumerable. Furthermore, the information contained in the shapes of curves “matters a great deal” (Horváth and Kokoszka,, 2012) and can be quite different from the information contained in their magnitudes. Therefore, it is important to develop methods that can isolate shape information and identify shape outliers from the rest of the data.

Many methods for identifying shape outliers rely on the notion of functional data depth. Functional data depth is a family of methods used to define centrality and induce a center-outward ordering on the sample functions (Liu et al.,, 1999; Zuo and Serfling,, 2000; Mosler and Polyakova,, 2012). All depth measures rank functions from most central (higher depth values) to least central (lower depth values) and are typically scaled to take on values in $[0,1]$.

Many outlier detection methods based on depth decompose total depth (or outlyingness) into magnitude and shape depth (or outlyingness). Arribas-Gil and Romo, (2014) proposed the Outliergram, a visualization tool for the shape and magnitude components of trajectories, using the half-space depth (Tukey,, 1977) and the band depth (Lopez-Pintado and Romo,, 2009). Later, Huang and Sun, (2019) introduced the total variation depth, which they decomposed into magnitude and shape components. In Dai and Genton, (2019), the authors pointed out that these methods, which rely on integrated depth, do not efficiently represent the centrality of functions. To remedy this, Rousseeuw et al., (2018) and Dai and Genton, (2019) simultaneously proposed the concept of directional outlyingess. Directional outlyingness has since been used as the basis of the functional outlier map (Rousseeuw et al.,, 2018) and the magnitude-shape plot tool (Dai and Genton,, 2018). Both of these methods again decompose their depth measures into magnitude and shape components, for the separate identification of magnitude and shape anomalies.

Other outlier detection methods proposed in the recent literature account for the geometry of the functions. Kuhnt and Rehage, (2016) developed the functional tangential angle (FUNTA) pseudo-depth based on the tangential angles of the intersections of the centered data. Nagy et al., (2017) proposed two modifications of previous depth notions to better identify shape outliers by emulating derivatives with multidimensional projections. Xie et al., (2017) separated the variability of functional data into amplitude and phase components, using the registration methods of Srivastava et al., (2011), and displayed this variability using independent boxplots for each component. Xie et al., (2017) showed how treating the phase and amplitude components of trajectories separately could greatly improve the detection rate of shape outliers. Their method, however, falls short of fully characterizing the shape distribution and instead relies on an optimization procedure to partially approximate the boundary of the inlier distribution.

Finally, there are methods based on Functional Principal Component Analysis (FPCA), which extract the features of normal shapes and detect anomalies by finding functions with abnormal features. These methods include a Step-wise Functional Outliers Detection test (Yu et al.,, 2012) and an FPCA score based distance test (Ren et al.,, 2017). FPCA based feature extraction can be quite powerful for detecting many types of shape anomalies, but FPCA is also known to be deficient when temporal variability is present (Srivastava et al.,, 2011; Tucker et al.,, 2013). Furthermore, these methods are not designed to explicitly identify shape anomalies separately from magnitude anomalies. Because we assume that both temporal variability and magnitude variability are present and are nuisance properties that mask shape, the problems and data considered here are quite different from the ones considered by FPCA methods. For those reasons, these methods were not included in the comparisons in Section 5.2.

In this paper, we introduce a new family of depth measures, called the elastic depths, based on the elastic shape distances used in Xie et al., (2017). We first use the elastic distances to directly measure the distance between the shapes of individual trajectories. This allows us to define a notion of data depth that appropriately captures the distribution of the trajectories’ shapes and allows for the explicit identification of shape outliers. Current literature falls short of fully representing the shape distribution and instead uses either a surrogate for shape or only approximates the shape distribution. Explicitly using the shape distribution, as we do in this paper, also confers several practical benefits, namely

Finding the inner quartile regions and outlier bound is trivial due to our depth based representation. The outlier bound does not require any optimization procedures and has only a single input parameter. We demonstrate in Section A.5 of the online supplement that the performance of our method is relatively insensitive to the value of the input parameter. 2. 2.

In Sections A.3 and 5.2, we empirically demonstrate the highly competitive and often superior outlier detection skill of the elastic depth based boxplots. We show that the elastic depths are consistently the top performing detection method across many different shape classes and that they can identify shape outliers in the presence of substantial translation and phase noise. 3. 3.

Shape distributions can easily be generalized to manifold valued trajectories, such as the hurricane trajectories (Figure 1), because the elastic depths are based purely on distance. Different manifolds merely require different distance metrics. In Section 6, we show the elastic depths applied to the Hurricane trajectories and in Section A.7 of the online supplement we show additional examples of the elastic depths applied to $\mathbb{R}$ and $\mathbb{R}^{2}$ valued data.

2 Background

2.1 Data Depth

Data depth is a general notion of measuring the centrality of observations with respect to a distribution. In the FDA literature, data depth is the dominant method used to define centrality and induce ordering on a function space. Given a distribution $P$ on a function space $F$, a depth function maps each trajectory $f\in F$ to a value in $[0,1]$, such that the closer a trajectory is to the center of $P$, the higher its depth value is. If this mapping is monotonic, that is, higher depth values necessarily mean higher centrality, then the depth function induces a center-outward ordering on the function space $F$ with respect to $P$. This makes depth a natural framework for evaluating the outlyingess of observations. High depth values mean an observation is very close to the center of $P$ so, conversely, low depth values mean an observation is very far from the center of $P$. Therefore, trajectories with extremely low depth values, as compared to the rest of the distribution, are likely to be outlying or anomalous.

2.2 Elastic shape analysis

Elastic shape analysis (ESA) is a collection of techniques for registering functional data through a phase-amplitude separation procedure and for performing statistical analysis on the separated phase and amplitude components (Srivastava et al.,, 2011; Kurtek et al.,, 2011; Tucker et al.,, 2013). Phase and amplitude represent two orthogonal components of a function’s variability. The amplitude component represents variability in shape, where shape refers to the properties of a function that remain unchanged under the shape preserving transformations: rotation, translation, scaling, and phase (Srivastava and Klassen,, 2016). The phase component represents the “domain” or “timing” variability of the trajectories. Because amplitude is invariant to these phase transformations, amplitude is distinct from the usual concept of magnitude. Magnitude measures the size of the observed realization of a trajectory while amplitude measures the size of the trajectories shape.

The distinguishing feature of ESA is the use of the Square Root Slope Function (SRSF) for registration (Kurtek et al.,, 2011). For real valued trajectories, the SRSF bijectively maps, up to an additive constant, a real valued function $f$ to its normalized gradient $f^{\prime}/\sqrt{|f^{\prime}|}$. Under ESA, two real valued trajectories are registered by elastically deforming the domain of one function such that the $L^{2}$ distance between the SRSFs of the two functions is minimized (Section 2.3). The amount of elastic deformation needed to register two functions is measured by the phase distance (Section A.1), while the residual $L^{2}$ distance between the SRSFs, post registration, defines the amplitude distance between them (Section 2.3). Together they are known as the elastic distances. The key insight of ESA is that by registering SRSFs, instead of trajectories directly, the amplitude distances are proper metrics and they are invariant to the shape preserving transformations. Thus, amplitude distance can be used to define the distance between the shapes of functions.

Later the Square Root Slope Velocity Function (SRVF) (Srivastava et al.,, 2011) was introduced to register $\mathbbm{R}^{n}$ valued trajectories and the Transported Square Root Slope Velocity Function (TSRVF) (Su et al.,, 2014) was introduced to register Riemannian manifold valued trajectories, such as trajectories observed on the unit sphere $\mathbbm{S}^{2}$. These notions allow us to calculate amplitude distances between multivariate functions and manifold valued functions respectively. We present the details for computing amplitude distances for $\mathbb{R}$ valued trajectories in Section 2.3. The details for $\mathbbm{R}^{n}$ valued and $\mathbbm{S}^{2}$ valued trajectories are deferred to the Appendix (Section A.1).

The advantages of the ESA approach to shape analysis have previously been shown in the works of Srivastava et al., (2011); Kurtek et al., (2011); Tucker et al., (2013); Su et al., (2014). ESA rigorously defines the shape space for a given class of trajectories and then defines a way to construct a proper distance metric on that shape space. The ESA based metrics are preserved under the shape preserving transformations: translation, scale, rotation, and reparameterization (phase). This improves theoretically over alternative shape metrics, such as Huang and Sun, (2016); Dai and Genton, (2018, 2019), that do not guarantee invariance or equivariance to shape transformations. The ESA framework is also general enough to apply to data observed in $\mathbbm{R}$, $\mathbbm{R}^{n}$, $\mathbbm{S}^{2}$, and any Riemannian manifold $\mathcal{M}$ that has an intrinsic metric. This is important for our motivating example, Atlantic hurricane trajectories (Figure 1), which are observed on the surface of a sphere.

2.3 Amplitude distance for $\mathbbm{R}$ valued functions

Let $F_{R}=\{f:[0,1]\mapsto\mathbb{R},f\text{ differentiable}\}$ be the class of differentiable trajectories on $[0,1]$ mapping to $\mathbbm{R}$. The Square Root Slope Function was introduced in Srivastava et al., (2011) as the following transformation on trajectories $f\in F_{R}$:

Definition 2.1**.**

Let $f$ be a differentiable trajectory in $F_{R}$, the Square Root Slope Function (SRSF) of $f$ is

[TABLE]

As was shown in their paper, the SRSF is a bijective mapping, up to an additive constant, from the space $F_{R}$ to the space of square integrable functions $L^{2}$. This means that for two functions $f,g\in F_{R}$, the norm on $L^{2}$

[TABLE]

where $q_{f},q_{g}$ are $f$ and $g$’s associated SRSFs, is a proper distance between $f$ and $g$ themselves. This norm is particularly important for shape analysis because it is phase invariant (Srivastava and Klassen,, 2016). That is, for any phase function $\gamma\in\Gamma$

[TABLE]

where $f\circ\gamma(t)=f(\gamma(t))$, $\forall t\in[0,1]$. Technical descriptions of phase functions $\gamma$ and phase space $\Gamma$ are deferred until Section A.1, but $\gamma$ functions essentially acts to deform the domain $[0,1]$.

Phase invariance means that Equation 1 is measuring some quantity that is independent of the representation, or phase, with which two functions are observed. This is only true if $f$ and $g$ share a common phase representation, so in order to find the amplitude distance between two arbitrary $f,g\in F_{R}$ we need to first place them in phase with each other. That is, we need to find some $\gamma^{*}\in\Gamma$ such that

[TABLE]

so that $||q_{f}-q_{g\circ\gamma^{*}}||_{2}$ measures the difference in their amplitudes. This can be more directly stated by defining the amplitude distance between $f$ and $g$ as in Srivastava et al., (2011):

Definition 2.2**.**

(Amplitude distance) Let $f$ and $g$ be two trajectories in $F_{R}$, then the amplitude distance between $f$ and $g$ is

[TABLE]

where $q_{f}$ and $q_{g\circ\gamma}$ denote the SRSF’s of $f$ and $g\circ\gamma$ respectively.

3 Elastic Depth

3.1 Definition of Elastic Depth

The exact analytic form of the elastic depths will greatly depend on the manifold on which the functional objects live. This is because the elastic depths are inherently going to be a robust summary of the distances between functional objects and the definition of distance between functional objects will inherently depend on the manifold on which they are observed. For instance, the distance between trajectories in $\mathbbm{R}^{n}$ is very different from the distance between trajectories on $\mathbbm{S}^{2}$. For the elastic distances to exist, however, it is only required that the data live on a Riemannian manifold (Section 3.1), such as $\mathbbm{R}^{n}$ or $\mathbbm{S}^{2}$, because the TSRVF of Su et al., (2014) can always be used to construct appropriate phase and amplitude distances. Therefore, we will only assume that our data live on a Riemannian manifold $M$ with an intrinsic metric, so the space of functions we consider is defined as

[TABLE]

The amplitude distance between two functions $f_{1},f_{2}\in F_{M}$ will generically be denoted as $d_{a}(f_{1},f_{2})$ and the phase distance between them as $d_{p}(f_{1},f_{2})$. The exact form of these distances is left unspecified because the amplitude distance is highly dependent on the manifold $M$. See Section 2.3 for the definition of amplitude distance for $\mathbb{R}$ valued trajectories and A.1 for the definitions of phase distance and amplitude distance for $\mathbbm{R}^{n}$ valued and $\mathbbm{S}^{2}$ valued trajectories. We now define the elastic depths for data observed on a manifold $M$ using the associated amplitude and phase distances. Let $P$ denote a distribution supported on the space $F_{M}$ and suppose we observe a function $f\in F_{M}$. We first introduce the idea of outlyingness, which describes the degree to which $f$ is an outlier relative to $P$. We further divide this concept into amplitude and phase outlyingness, using the amplitude and phase distances respectively. This is done to separately quantify the shape outlyingness and phase outlyingness of $f$ relative to $P$. Amplitude and phase outlyingness are respectively denoted as $O_{a}$ and $O_{p}$ and are defined as

[TABLE]

where $X$ is a random function in $F_{M}$ and drawn from the distribution $P$. The outlyingness functions $O_{a}$ and $O_{p}$ robustly summarize the pairwise distances between $f$ and all other functions $X\in F_{M}$. These two functions define a measure of outlyingness such that if $O_{a}(f,P)$ is large then $f$ is generally dissimilar in amplitude from other functions $X\in F_{M}$, with respect to the distribution $P$. Likewise, if $O_{a}(f,P)$ is small then $f$ is similar in amplitude to other functions $X\in F_{M}$.

To convert $O_{a}$ and $O_{p}$ into depth functions we invert them with the type B depth construction of Zuo and Serfling, (2000):

[TABLE]

$D_{a}(f,P)$ and $D_{p}(f,P)$ are respectively called the amplitude depth and phase depth of $f$ with respect to $P$. Together we denote them the elastic depths. The purpose of inverting the outlyingness functions in this manner is to create bounded measures of centrality, i.e. depths, on the amplitude and phase spaces associated with $F_{M}$. When depths, such as $D_{a}$ and $D_{p}$, satisfy the properties outlined in Section 3.2, they provide a non-parametric and moment free characterization of the distribution $P$. Larger depth values indicate higher centrality and low outlyingness, while lower values indicate higher outlyingness. Thus, $D_{a}$ and $D_{p}$ provide a simple and rigorous way to identify outliers based on the underlying distribution $P$.

3.2 Properties

Within the depth literature there have been many desirable properties discussed for both multivariate and functional data depths; see Zuo and Serfling, (2000) and Mosler and Polyakova, (2012) for comprehensive reviews. These properties ensure that a depth function properly measures the notion of depth or centrality. For instance, a depth function needs to be location and scale invariant (or equivariant) and it should decrease monotonically from a natural point of symmetry. Since our depth is purely for functional data we concentrate on the central properties of Mosler and Polyakova, (2012). These properties are established for the amplitude depths because amplitude is the primary concern of shape analysis.

The elastic depths are based on proper distance metrics so they inherit certain properties such as translation invariance and scale equivariance automatically. On some manifolds, such as $\mathbbm{R}^{2}$, scale equivariance can be promoted to scale invariance because the trajectories are constrained to live on an $L^{2}$ ball. Invariance to simultaneous reparameterization (simultaneous phase invariance) was shown in Srivastava et al., (2011) for amplitude distances between $\mathbbm{R}$ and $\mathbbm{R}^{n}$ valued trajectories and then later extended to $\mathbb{S}^{2}$ valued trajectories in Su et al., (2014). Consequently, the amplitude depths are also invariant to simultaneous reparameterization.

Other properties, such as phase invariance, maximality of the center, and convex level sets are essential for shape anomaly detection but are not simple corollaries of the amplitude distance. We outline these properties, as they apply to amplitude depth, below. All proofs are deferred to the online supplement Section A.2.

Proposition 3.1** (Phase invariance).**

Let $\Gamma$ be the space of warping, or phase, functions defined in Section A.1 and let $\gamma\in\Gamma$. Let $F_{M}$ be the space $M$-valued differentiable functions as in Section 3.1, let $f\in F_{M}$ and suppose $P$ is a distribution supported on $F_{M}$. Then

[TABLE]

where $D_{a}(\cdot,P)$ is the amplitude depth of trajectories on $F_{M}$ with respect to $P$.

This property is unique to the elastic depths and ensures that the amplitude depths are invariant to the phase under which each trajectory is observed. This property, in conjunction with translation and scale invariance (equivariance), means that the amplitude depth is invariant to the shape preserving transformations. We can, therefore, say that amplitude depths are appropriately capturing our definition of shape.

Proposition 3.2** (Maximality of the center).**

Let $F_{M}$ be the space $M$-valued differentiable functions as in Section 3.1, let $f\in F_{M}$ and suppose $P$ is a distribution supported on $F_{M}$. A trajectory $s\in F_{M}$ is the amplitude Fréchet median of $P$ if and only if $s=\arg\max_{f\in F_{M}}D_{a}(f,P)$, where $D_{a}(\cdot,P)$ is the amplitude depth of trajectories on $F_{M}$ with respect to $P$.

Maximality of the center guarantees that the maximizer of the amplitude depths, denoted the amplitude depth median, is the actual Fréchet median of the distribution. The Fréchet median is the trajectory that minimizes the median distance between itself and all other points in the space. This property ensures that the amplitude depths start their ordering from the true amplitude center of the distribution.

Proposition 3.3** (Convex level sets).**

Let $F_{M}$ be as in Section 3.1, let $f\in F_{M}$ and suppose $P$ is a distribution supported on $F_{M}$. Let $D_{a,\alpha}(P)=\{f\in F_{M}:D_{a}(f,P)\geq\alpha\}$ be the upper level sets for the amplitude elastic depth for all $\alpha\in[0,1]$. Then $D_{a,\alpha}(P)$ is a convex set. Similarly the upper level sets for the phase elastic depth $D_{p,\alpha}(P)=\{f\in F_{M}:D_{p}(f,P)\geq\alpha\}$ are convex for all $\alpha\in[0,1]$.

Convexity of the level sets implies that depths decrease monotonically from the center of the distribution. In conjunction with Maximality of the Center, level set convexity guarantees that the elastic depths are measuring centrality in amplitude space and phase space. This property further distinguishes the elastic depths from previous depth notions because they do not directly characterize centrality in the appropriate shape spaces. We use these convex level sets as the theoretical basis for the construction of the depth boxplots (Section 4) and for depth thresholding (Section 4.2).

3.3 Estimating Elastic Depths

As in Section 3.2, let $F_{M}$ be the space of differentiable functions on the Riemannian manifold $M$ and let $P$ represent a distribution supported on $F_{M}$. Suppose we observe $f_{1},...,f_{n}\sim P$. The amplitude and phase depths of each $f_{i}$, $i\in 1,...n$, can be estimated empirically using their respective sample outlyingness functions. The sample amplitude and phase outlyingness functions are respectively denoted as $O_{a,n}$ and $O_{p,n}$ and are defined as:

[TABLE]

where $P_{n}$ denotes the empirical distribution of the functions $f_{1},...,f_{n}$. Using the same construction as before, we invert the sample outlyingness functions into sample depths

[TABLE]

for amplitude and phase respectively. The following proposition asserts the uniform consistency of this depth estimator.

Proposition 3.4** (Uniform Consistency).**

Let $F_{M}$ be as in Section 3.1, suppose $P$ is a distribution supported on $F_{M}$, let $f_{1},...,f_{n}\sim P$, and let $P_{n}$ represent the empirical distribution of the sample. Then

[TABLE]

where $D_{a}(\cdot,P)$ is the amplitude depth of trajectories on $F_{M}$ with respect to $P$, $D_{a,n}(f,P_{n})$ is the amplitude depth’s empirical counterpart, $D_{p}(\cdot,P)$ is the phase depth of trajectories on $F_{M}$ with respect to $P$, and $D_{p,n}(f,P_{n})$ is the phase depth’s empirical counterpart.

4 Identifying outliers

Data depth is a natural framework for outlier detection because it provides a center-outwards ordering of the data. Functions with very low depth values are strong candidates for outliers because they are statistically far from the center of the distribution. As mentioned in Section 1, there have been many methods, many based on functional depth in some way, for detecting shape anomalies proposed in the literature. These methods typically construct an outlier cutoff boundary on either the depths or the functions and classify any trajectory as an outlier if it exceeds these bounds. In the next two sections, we introduce two simple ways of defining an outlier cutoff point based on elastic depth.

4.1 Depth Boxplots

The first method we introduce is called the Depth Boxplot, which is a half-boxplot constructed on the elastic depths directly. We showed in Section 3.2 that the elastic depths decrease monotonically from their unique center, as trajectories become more outlying, so using the depths directly does not incur a loss of outlyingness information. Additionally, unlike methods that place bounds on the observed data, using a boxplot on the depth values circumvents the problem of shape outliers being masked due to scale, translation, and phase variability. This is because the boundaries of a depth boxplot correspond to entire central regions on the shape space of functions, and not merely central regions on the projections of functions onto $\mathbbm{R}$ (or $\mathbbm{R}^{2}$ or $\mathbbm{S}^{2}$).

Algorithm 1 describes how to construct the amplitude depth boxplot and how amplitude anomalies are identified with the whisker $c$. Phase anomalies can similarly be defined by substituting amplitude depths for phase depths.

The boxplot created in Algorithm 1 consists of the following three pieces: The median, the IQR, and the whisker $c$ (Figure 2). The median of the boxplot is the largest depth, because as was shown in Section 3.2, the largest depth corresponds to the median of the distribution. The IQR is the 50% central region because, as per the IQR of univariate data, this range contains the inner 50% of the data. Most importantly is the whisker value $c$, which determines which trajectories are considered outliers. Any trajectory with an amplitude depth of less than $c$ is considered an anomaly because it is statistically too far from the rest of the data.

The whisker $c$ is determined by a multiplier or inflation factor $k$. The quantity $k$ is a free parameter that must be set to detect anomalies. In classical univariate boxplots, $k=1.5$, so as to achieve approximately 99.3% coverage of the boxplot on Gaussian data. This guarantee does not necessarily extend to functional data, but we find empirically that $k\in[1.5,2]$ works well to separate outliers from inliers as long as the sampling frequency is high enough to fully represent the functional data. We investigate the depth boxplot’s dependency on $k$ numerically in Section A.5 of the online supplement and find that detection performance is fairly robust to $k$ for a wide range of values and across many types of data. More details on the choice of $k$ can also be found in Section 5.2.

4.2 Depth Thresholding

Boxplots and other hard cutoffs are not the only way to investigate shape anomalies. The quantiles of the depth distribution itself can also be used to investigate the most extreme data. The elastic depths induce a proper center-outward ordering of the trajectories, so the most extreme, i.e. smallest, depth values correspond to the most extreme trajectories. Therefore, if we wanted to view the 5% most extreme functions, we could simply select the functions with the 5% smallest elastic depth values. More generally, the $100\times p$% most outlying functions have depth values below the $(1-p)^{th}$ quantile of the depth’s distribution. This type of thresholding is quite useful in exploratory analysis for comparing, say, the 1%, 5%, and 10% most outlying shapes with the 1%, 5%, and 10% most inlying shapes.

The limitation of thresholding is that it will always select $p$% of the data to be outlying, so as an anomaly detector it is insufficient on its own. However, it can be paired with the depth boxplots to produce a more robust depth boxplot. Algorithm 2 extends algorithm 1 to include a thresholding parameter $p$ so that to be considered an outlier, a function must have a depth value below the whisker $c$ and below the $(1-p)^{th}$ quantile of the depth’s distribution.

The purpose of $p$ is mainly to control the number of false positives when detecting outliers. If $p=0.95$, then at most 5% of the trajectories will be considered outlying, no matter what the whisker value is. While the whisker is generally sufficient for achieving good coverage of the depth distribution and detecting anomalies, there are situations, such as low sampling frequency (Section A.5 of the online supplement), where the whisker can fall short unless the multiplier $k$ is set higher. In these situations, $p$ will act to effectively increase the $k$ so that better coverage is achieved.

5 Simulation Study

A simulation study was conducted to comprehensively assess the performance of the elastic depths and the depth boxplots. We compared our method against nine other shape anomaly detectors: the Outliergram (OG) (Arribas-Gil and Romo,, 2014), Sequential Transformations (ST-T1, ST-T2, ST-D1) (Dai et al.,, 2018), the Functional Outlier Map (FOM) (Rousseeuw et al.,, 2018), Total Variation Depth (TVD) (Huang and Sun,, 2016), the Magnitude-Shape (MS) plot (Dai and Genton,, 2018), the Robust Functional Tangential Angle Pseudo-depth (rFUNTA) (Kuhnt and Rehage,, 2016), Order Extended Integrated Depth (FDJ and IDJ) (Nagy et al.,, 2017), Geometric boxplots (GEOM) (Xie et al.,, 2017), and Directional Outlyingness (DIR) (Dai and Genton,, 2019). In the following sections comparisons to TVD, MS, DIR, and GEOM are included while the rest are deferred to the appendix. These four methods were consistently the strongest competitors across each of the outlier models. We describe their implementation here briefly.

The TVD outliers were found using the detectOutlier function in the TVD R package with an empFactor = 1.5. MS outliers were found by computing the MO and VO quantities then using the cerioli2010.irmcd.test function in the CerioliOutlierDetection R package to compute the boundary with a coverage probability of 99.3%. DIR outliers were found using the authors dir.out function with Mahalanobis distance and the default parameters fac $=0.154$, and cutoff $=6.91$. GEOM outliers were found using the AmplitudeBoxplot function in the fdasrvf R package using $k=1$. Implementation details for the detectors that are deferred to the appendix have likewise been deferred to the appendix (Section A.4 of the online supplement). Elastic depths (ED) outliers were identified using the depth boxplots with $k=1.8$. Boxplots were computed on both the amplitude and phase depths separately. The results using amplitude depth are denoted as ED-A and the results using phase depths are ED-P.

5.1 Simulation Design

We define seven different shape outlyingness scenarios to test the effectiveness of the above shape outlier detectors. Each of these scenarios is represented by one of the seven models detailed below. The first six correspond to amplitude (shape) outliers while the seventh is for phase outliers.

Model 1 (Amplitude Increase): Main model: $X(t)=\sin(5\pi t)+4t+e(t)+\delta$ and Contamination model: $X(t)=4\sin(5\pi t)+4t+e(t)+\delta$, where $t\in[0,1]$, $e(t)$ is a centered Gaussian process with covariance function $\gamma(x,x^{\prime})=\exp\{-(x-x^{\prime})^{2}/0.5\}$, and $\delta\sim N(0,1)$ is a random additive translation term. The purpose of $\delta$ is to shift each curve by a random amount so as to mask shape outliers that could accidentally be identified as magnitude outliers. 2. 2.

Model 2 (Amplitude Decrease): Main model: $X(t)=\sin(5\pi t)+4t+e(t)+\delta$ and Contamination model: $X(t)=\frac{1}{6}\sin(5\pi t)+4t+e(t)+\delta$, where $t\in[0,1]$, and $e(t)$ is the Gaussian process from Model 1. 3. 3.

Model 3 (Mixed Polynomials): Main model: $X(t)=t^{3}-2t^{2}+0.5t+e(t)$ and Contamination model: $X(t)=2t^{3}+t^{2}-0.5t+e(t)$. 4. 4.

Model 4 (Covariance change): Main model: $X(t)=\sin(5\pi t)+4t+e_{1}(t)+\delta$ and Contamination model: $X(t)=\sin(5\pi t)+4t+e_{2}(t)+\delta$, where $t\in[0,1]$ and $e_{1}(t)$ and $e_{2}(t)$ are centered Gaussian processes with covariance functions $\gamma(x,x^{\prime})=\exp\{-(x-x^{\prime})^{2}/50\}$ and $\gamma(x,x^{\prime})=\exp\{-(x-x^{\prime})^{2}/2\}$, respectively. 5. 5.

Model 5 (Frequency Increase): Main model: $X(t)=\sin(2\pi t)+4t+e_{1}(t)+\delta$ and Contamination model: $X(t)=\sin(12\pi t)+4t+e(t)+\delta$ where $t\in[0,1]$ and $e(t)$ is the Gaussian process from Model 1. 6. 6.

Model 6 (Jump contamination): Main model: $X(t)=\sin(5\pi t)+4t+e(t)+\delta$ and Contamination model: $X(t)=\sin(5\pi t)-2\mathbbm{1}_{(t<T)}+3\mathbbm{1}_{(T\leq t)}+4t+e(t)+\delta$, where $t\in[0,1]$ and T is distributed uniformly on $[0.4,0.6]$. 7. 7.

Model 7 (Phase Contamination): Main model: $X(t)=\sin(5\pi t)+4t+e(t)+\delta$ and Contamination model: $X(t)=\sin(5\pi\gamma(t))+4\gamma(t)+e(\gamma(t))+\delta$, where $t\in[0,1]$ and $\gamma$ is a random phase function from $\Gamma$. The functions $\gamma$ are generated from the first two Fourier basis functions with random amplitudes distributed as $N(0,\sigma)$ on the tangent space to the unit Hilbert sphere. We use $\sigma=6$ to impose a large amount of phase variability on the contamination model (Figure 5).

Each of these models, except model 7, was then further contaminated with two additional sources of noise: compositional (phase) noise and magnitude outliers. Compositional noise was added by composing each trajectory with a random phase function generated by the rgam function in the fdasrvf package with sigma = 0.1. Magnitude outliers were added by randomly shifting $10\%$ of the generated functions by $\pm 10$. Model 7 was only contaminated with magnitude outliers because adding phase noise would destroy the difference in phases that we are trying to detect. These two noise sources introduce a level of realism to our simulations because nuisance phase and magnitude outlyingness are often present when analyzing shapes. The base amplitude outlier models without additional phase and magnitude noise are pictured in Figure 3.

5.2 Contamination by Multiple Anomalies

We considered the case when 10% of the data is outlying in shape. We compared the performance of the detection methods on the seven outlier models using the $F_{1}$ score (Chinchor,, 1992) for outlier classification. The $F_{1}$ score is a comprehensive measure of classification accuracy that considers both the precision (positive predictive value) and the recall (true positive rate) of a detection method. A method that perfectly classifies all outliers as outliers and all inliers as inliers will have an $F_{1}$ score of 1. Methods that do not perfectly classify will have $F_{1}$ scores less than 1. The $F_{1}$ score was traditionally defined as the harmonic mean of precision and recall, but it can also be expressed in terms of the more familiar True Positive (TP), False Negative (FN), and False Positives (FP) quantities:

[TABLE]

90 inlying trajectories and 10 outlying trajectories were sampled from the main model and contamination model respectively. Compositional noise and magnitude outliers were again added to each of the models, except for model 7 (Phase Contamination) where only magnitude outliers were added. Trajectories were sampled on an equidistant 30 point grid over $[0,1]$ and 1000 simulations were performed for each of the six models. The results for the top models are summarized in Figure 4. Full results for all considered models are available in the appendix.

Figure 4 shows that, across the six amplitude outlier models, the elastic depths maintained the highest average $F_{1}$ score. The amplitude depth based boxplots have an average $F_{1}$ score of around 0.95 - 1.0, indicating that in each scenario, regardless of outlier or inlier type, the amplitude depths can achieve near-perfect detection. Consistently high performance of the elastic depths is notable because existing methods, while strong in some cases, suffer major losses of power in others. For instance, on some models, such as model 3 (mixed polynomials) and model 4 (covariance change), the amplitude depth based boxplot was the only method able to detect the shape outliers consistently. Even GEOM, which uses the elastic distances, was unable to consistently separate these outliers from the inliers. Together, the results show that because the elastic depth based boxplots use the shape distribution, albeit indirectly via data depth, they can consistently and skillfully detect a wide variety of shape outliers. They do not generally suffer a loss of power due to compositional (phase) noise, translation noise, presence of magnitude outliers, or even inlier and outlier type.

Figure 5 shows that, across the top phase outlier detection models, ED-P, the phase counterpart of ED-A, has near-perfect detection skill. The TVD, however, maintains a slightly higher $F_{1}$ score. Our conjecture is that the data generation process may have incidentally induced non-phase based differences into the phase outliers. TVD, MS, and DIR could potentially take advantage of this non-phase information to improve their score, whereas ED-P and GEOM could not because they only use phase information. See a small simulation example in Section A.6 in the online supplement.

We also investigated the sensitivity of the boxplots to the parameter $k$ via a simulation study in Section A.5 of the online supplement. Overall, we found that the coverage of the boxplot on the inlying uncontaminated data is insensitive to the value of $k$. Across each model, the coverage of the boxplots steadily increases from about 95% at $k=1$ to around 100% at $k=3$. Some models required higher values of $k$ to achieve the desired 99% coverage, which we found to be due to under sampling, i.e. sampling below the Nyquist rate, of the observed trajectories. We recommend avoiding this potential issue by ensuring adequate sampling of the trajectories when possible or using $k\approx 2$ when this is not possible.

6 Hurricane Trajectories

Our motivating example (Figure 1) comes from the National Hurricane Center’s (NHC) Atlantic Hurricane Database (HURDAT2) (Landsea and Franklin,, 2013). The NHC assimilates all observations, real time and post-storm, for each tropical cyclone to estimate and record its characteristics and path across the Atlantic Ocean. The HURDAT2 database contains records for 979 tropical cyclone paths of various lengths, shapes, sizes, orientations, and placements. We only consider storms with at least 25 observations because only those paths had sufficient time to develop. The storms were then further subset to include only those originating in the ocean and north of South America.

The typical path of a hurricane is “U” shaped, starting in Africa then cutting across the eastern United States and finally heading back east towards Europe. Due to the vast distances hurricanes travel it would be inappropriate to treat them as lying on a Euclidean plane. Instead, we consider them as trajectories on the surface of a unit sphere $\mathbb{S}^{2}$. We used the elastic depth boxplots with a depth threshold of 0.05 to limit the number of amplitude outliers to fourteen; the top four of which are pictured in Figure 6.

Each of the top four outliers is markedly different from the standard “U” shape. They exhibit an atypical spiraling behavior as they meander across the Atlantic. The identification of shape outliers helps climate scientists further investigate what causes the trajectories to be anomalous. It can be very important for improving the accuracy of hurricane prediction algorithms if the dynamics which produce anomalies are well understood. Further data examples on $\mathbbm{R}$ and $\mathbbm{R}2$ valued trajectories can be found in section A.7 of the online supplement.

7 Discussion

In this paper, we proposed a new class of functional depths based on the elastic distance metrics and showed how they may be used to detect shape outliers. The theoretical properties of our new elastic depth were investigated, and it was shown that they satisfy most key properties required of a depth metric. These include translation, scale, rearrangement, phase invariance (equivariance), maximality of the center, convex level sets and monotonicity from the center. Rearrangement invariance and phase invariance were particularly crucial for detecting shape anomalies because they allowed the amplitude depth to measure centrality independent of phase.

We demonstrated the empirical performance of our method together with nine competing methods using extensive simulation studies. It was shown that our method attains the overall highest average $F_{1}$ score across all models. On each of the six amplitude models the amplitude depth based boxplots had an average $F_{1}$ score of around 0.95 - 1.0, indicating that in each scenario, regardless of outlier or inlier type, the amplitude depths achieved near-perfect detection on average. On some models, such as model 3 (mixed polynomials) and model 4 (covariance change), the amplitude depth based boxplot was the only method able to detect the shape outliers consistently. Together, these results demonstrate the power of using the shape distribution, albeit indirectly via data depth, to detect shape outliers.

The simulation results, of course, depend on the boxplot multiplier $k$. We recommend setting $k=2$ as the default value. We found through empirical studies (See online supplement Section A.5) that any $k$ value between 1.5 and 2.25 provides nearly the same level of detection skill. The value $k=2$ had the most favorable trade-off between false positives and false negatives across all of the outlier models in the simulation. Larger values (up to around 2.5-2.8) are helpful when the sampling frequency of the trajectories is too low (See online supplement Section A.5). If the sample size is large then the depth quantiles, $p$, can also be used to help set a minimum value for $k$. That is, choose $k$ large enough to make the whisker value less than the $(1-p)$’th quantile of the amplitude depth. This guarantees that the shapes of the inner $p\times 100$% of trajectories will be considered as inliers.

Finally, we showed how the elastic depths may be used to identify shape outliers in functional data observed on the unit sphere $\mathbb{S}^{2}$. We used the HURDAT2 hurricane track database and identified the four most shape outlying trajectories. We found that these trajectories’ paths were remarkably different from the standard U-shaped paths that hurricanes normally follow. Further applications in the online supplement A.7 demonstrate the elastic depths on $\mathbbm{R}$, $\mathbbm{R}^{2}$ valued functional data. These data examples illustrated the simplicity and consistency with which the elastic depths may be applied, regardless of the underlying geometry of the space.

Acknowledgments

This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government. This work was supported by the Laboratory Directed Research and Development program at Sandia National Laboratories, a multi-mission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525.

Supplementary Materials

Supplementary results:

PDF file “Supplementary results” containing detailed descriptions of the elastic distances, proofs for all results in Section 3.2, additional simulations, and two additional real data examples. (pdf file).

Source code:

Zip file “Source code” containing code to produce all figures and results in the manuscript and supplementary file. (zip file)

R-package for Elastic Depth:

R-package “elasticdepth” containing code to compute the elastic depths and depth boxplots for $\mathbb{R}$, $\mathbb{R}^{2}$, and $\mathbb{S}^{2}$ valued trajectories. Available through GitHub: trevor-harris/elasticdepth. (GNU zipped tar file)

Appendix A Appendix

A.1 Amplitude distance

In this section, we will introduce the explicit definitions and mathematical motivation for the amplitude distances mentioned in Section 2.2. Later parts of this section rely on concepts from topology and differential geometry such as tangent spaces and parallel transport. A thorough review of these topics, as applied to the problem of shape analysis, can be found in Chapter 3 of Srivastava and Klassen, (2016). Here we will only provide conceptual introductions as they are needed to define amplitude metrics.

Phase Space

To define the concept of amplitude we first need to define the concept of phase, because amplitude is the properties of a function that are independent of phase. The phase space of the unit interval $[0,1]$ is defined as

[TABLE]

This diffeomorphic constraint gives rise to the notion of elasticity because the elements of $\Gamma$, i.e. phase functions, can only smoothly stretch and contract portions of the unit interval so that it maps back to itself. Phase is generally thought of as the representation of a trajectory because any trajectory with domain $[0,1]$ can be warped by a phase function to appear differently. The amplitude will be taken to be those features of a trajectory that remain unchanged under any possible warping. In the following sections, we will consider functions with domains on $[0,1]$ which take values in $\mathbb{R}^{2}$ and the unit sphere $\mathbb{S}^{2}$ and develop proper distance metrics on their amplitudes.

Amplitude distance for $\mathbbm{R}^{n}$ valued functions

Let $F_{\mathbbm{R}^{n}}=\{f:[0,1]\mapsto\mathbb{R}^{n},f\text{ differentiable}\}$ be the class of differentiable trajectories on $[0,1]$ mapping to $\mathbbm{R}^{n}$ with $n\geq 2$. In higher dimensional Euclidean space ($n\geq 2)$, the scale, rotation, and phase of a trajectory have to be accounted for in order to isolate its shape. Scale variability is removed by standardizing each trajectory to have length one, that is each trajectory is divided by the $L^{2}$ norm of its gradient:

[TABLE]

where $||f^{\prime}||=\sqrt{\int_{0}^{1}f^{\prime}(t)^{2}dt}$. This operation can intuitively be thought of as standardizing the “zoom level” or “magnification” of each trajectory. In Section 2.3, the $\mathbbm{R}$ valued trajectories were not length standardized because, on $\mathbbm{R}$, scale and amplitude are essentially the same after trajectories are aligned.

Accounting for rotation variability requires introducing the space of rotation matrices, the Special Orthogonal Group $SO(n)$, and another level of optimization. The space $SO(n)$ is defined as the group of orthogonal matrices with determinant one. For trajectories in $\mathbbm{R}^{2}$, we use $SO(2)$ to define rotations around a point, in $\mathbbm{R}^{3}$ the group $SO(3)$ defines rotations around a line, etc. The action of $SO(n)$ on a trajectory $f$ is denoted as $O(f)$ and is defined pointwise as

[TABLE]

where $Of(t)$ represents standard matrix multiplication of the $n\times n$ matrix $O$ and the $n\times 1$ vector $f(t)$. See Vladimir, (1994) for more details on the Orthogonal groups and their properties. Finding the optimal rotation matrix in $SO(n)$ is done alongside finding the optimal phase function when computing the amplitude distance.

To represent phase variability in $\mathbbm{R}^{n}$ valued function we can again use elements of the space $\Gamma$. This is because $\Gamma$ is only defined with respect to the domain $[0,1]$ and not the range. As in the real valued function case, we want a distance that is invariant to simultaneous reparameterizations so we use the Square Root Velocity Function (SRVF), which was introduced in Srivastava et al., (2011) as the following transformation on trajectories $f\in F_{R}^{n}$:

Definition A.1**.**

Let $f$ be a differentiable trajectory in $F_{R}^{n}$, the Square Root Slope Function (SRVF) of $f$ is

[TABLE]

where $|f^{\prime}(t)|$ is the absolute value of $f^{\prime}$ at $t$.

This definition is effectively the multivariate generalization of Definition 2.1. The only difference is that, due to the length restriction on $f$, this transformation maps trajectories onto the $L^{2}$ ball of radius one instead of $L^{2}$ space itself (Srivastava and Klassen,, 2016). Without this restriction the SRVF would map to $L^{2}$ and be an exact generalization of the SRSF.

This definition is nearly identical to that of 2.1, except that the trajectory $f$ takes values in $\mathbbm{R}^{n}$. Due to the length constraint on $f$, this transformation maps the length constrained functions onto the $L^{2}$ ball of radius one, instead of $L^{2}$ space itself (Srivastava and Klassen,, 2016). This means the norm on the SRVFS is instead the arc length distance on $L^{2}$ spheres which is defined as

[TABLE]

where $q_{f}=SRVF(f)$ and $q_{g}=SRVF(g)$ for two trajectories $f,g\in F_{R^{n}}$ and $\langle q_{f}(t),q_{g}(t)\rangle$ denotes the inner product of the vectors $q_{f}(t),q_{g}(t)$. To convert this into an amplitude distance we need to place $f$ and $g$ in phase and rotation with each other. This optimization is summarized as the following amplitude distance on $\mathbbm{R}^{n}$ valued trajectories (Srivastava et al.,, 2011):

Definition A.2**.**

Let $f$ and $g$ be two trajectories in $F_{R^{n}}$, then the amplitude distance between $f$ and $g$ is

[TABLE]

where $q_{f}$ and $q_{O(g\circ\gamma)}$ denote the SRVF’s of $f$ and $O(g\circ\gamma)$ respectively.

Amplitude distance for $\mathbb{S}^{2}$ valued functions

Let $F_{\mathbbm{R}^{2}}=\{f:[0,1]\mapsto\mathbb{S}^{2},f\text{ differentiable}\}$. For differentiable functions constrained to live on a sphere, rotation about a point and scaling are not possible. However, the direct comparison of function gradients (SRSVs) is also not possible without additional steps (Srivastava and Klassen,, 2016). The gradient of a function at a particular point lives in a vector space perpendicular to $\mathbbm{S}^{2}$ known as the tangent space. The tangent space represents all possible “directions” in which a trajectory can pass through a given point, and can be viewed as a flat plane touching $\mathbbm{S}^{2}$ only at that point.

On a nonlinear manifold, such as $\mathbbm{S}^{2}$ the tangent space changes depending on which point on the manifold it touches. This means the distance between gradients (SRVFs) is not well defined because they don’t live in the same space. To remedy this issue Su et al., (2014) introduced the Transported Square Root Velocity Function (TSRVF), which uses the idea of parallel transport from differential geometry to make the tangent spaces of two trajectories on $\mathbbm{S}^{2}$ comparable. Roughly speaking, parallel transport describes how to transform the tangent space of one location into the tangent space of another location, provided the two locations can be connected by a smooth path (Kobayashi and Nomizu,, 1996). On the unit sphere $\mathbbm{S}^{2}$ parallel transport is a straightforward analytic calculation (Srivastava and Klassen,, 2016).

Definition A.3**.**

Let $f$ be a differentiable trajectory in $F_{\mathbb{S}^{2}}$ and let $c\in\mathbb{S}^{2}$. Then $\forall t\in[0,1]$ the parallel transport of $f^{\prime}(t)$ to the tangent space of the point $c$, along the shortest available path, is defined as

[TABLE]

This definition means that $f^{\prime}(t)$, which is tangent to $f(t)$, can be converted into tangent vector at any point $c$ on the sphere. More importantly this means that the entire gradients of two trajectories $f,g$ can both be moved into the tangent space of a single point. Once both trajectories completely share a common tangent space, the amplitude distance between $f,g$ is computable. This leads to the notion of the TSRVF in (Su et al.,, 2014):

Definition A.4**.**

Let $f$ be a differentiable trajectory in $F_{M}$ for some Riemannian manifold $M$ with norm $||\cdot||$. The Transported Square Root Vector Field (TSRVF) of $f$ at the point $f(t)$ to the point $c\in M$ is

[TABLE]

Where $f^{\prime}(t)_{f(t)\mapsto c}$ represents the parallel transport of the tangent vector $f^{\prime}(t)$ from $f(t)$ to $c$.

To define an amplitude distance, a common point $c\in\mathbbm{S}^{2}$ must be chosen to transport each function to. As Su et al., (2014) showed, the choice of $c$ only minorly impacts the computed distances and does so in a consistent way. We also have to account for phase variability as we did for $\mathbbm{R}$ and $\mathbbm{R}^{n}$ valued trajectories. This leads to the following definition of amplitude distance for trajectories on $\mathbbm{S}^{2}$ (Su et al.,, 2014):

Definition A.5**.**

Let $f$ and $g$ be two trajectories in $F_{S^{2}}$, then the amplitude distance between $f$ and $g$ is

[TABLE]

where $h_{f}$ and $h_{g\circ\gamma}$ denote the TSRVF’s of $f$ and $g\circ\gamma$ respectively to some common point $c\in\mathbbm{S}^{2}$.

Phase distance

The phase distance is the same for each of the three manifolds and is, in fact, the same for all univariate functions mapping to any Riemannian manifold. This is because the phase space $\Gamma$ is not defined with respect to the range of the functions, $M$, but only to the domain $[0,1]$. To define phase distance we use the optimal $\gamma$ that defines the amplitude distance. The phase space $\Gamma$ is a nonlinear manifold with no known geometry so we use the SRSF to map $\Gamma$ to a known geometry. Phase functions are positive for all $t\in[0,1]$ and $||q_{\gamma}||=\int_{0}^{1}\sqrt{\gamma^{\prime}(t)}^{2}dt=1$, so the SRSF maps $\Gamma$ onto the positive orthant of a unit Hilbert Sphere. Thus the phase distance is defined as

[TABLE]

where $I(t)=t$ is the identity function and $\langle\sqrt{I^{\prime}(t)},\sqrt{\gamma^{\prime}(t)}\rangle$ is the inner product between the two vectors $\sqrt{I^{\prime}(t)}$ and $\sqrt{\gamma^{\prime}(t)}$ (Srivastava and Klassen,, 2016). The metric $d_{p}(f_{1},f_{2})$ is essentially measuring the amount of elastic deformation needed to compare the amplitudes of $f_{1}$ and $f_{2}$.

A.2 Proofs

Proof of Proposition 3.1 (Phase Invariance)

Proof.

Let $f\in F_{M}$ and let $\gamma\in\Gamma$. Then

[TABLE]

by phase invariance of the amplitude distances (Kurtek et al.,, 2011). ∎

Proof of Proposition 3.2 (Maximality of the Center)

Proof.

The Fréchet median $m$ of the amplitudes is defined as the trajectory which minimizes the median amplitude distance between itself and all other trajectories in $F_{M}$ (Srivastava and Klassen,, 2016). That is, $m=\arg\min_{f\in F_{M}}\text{median}\{d_{a}(f,X):X\in F_{M}\}$. Therefore if $s=\arg\max_{f\in F_{M}}D_{a}(f,P)$ then:

[TABLE]

Thus the maximizer of the amplitude depth is the Fréchet median of the amplitudes. ∎

Proof of Proposition 3.3 (Convex Level Sets)

Proof.

Let $f_{1}$ and $f_{2}$ be in $F_{M}$ with amplitudes $[q_{1}]$ and $[q_{2}]$ in $D_{a,\alpha}(P)$. Let $f\in F_{M}$ such that the amplitude of $f$ is $[q]=\lambda[q_{1}]+(1-\lambda)[q_{2}]$ for some $\lambda\in[0,1]$. The amplitude distance between $f$ and a random $X\sim P$ in $F_{M}$ with amplitude $[q_{X}]$ can be upper bounded as follows

[TABLE]

by the convexity of the amplitude distance. Since this is a convex combination of positive real numbers, $d(\cdot,\cdot)$, we have that

[TABLE]

Since $D_{a}(f_{1},P)\geq\alpha$ and $D_{a}(f_{2},P)\geq\alpha$, by virtue of $[q_{1}]$ and $[q_{2}]$ in $D_{a,\alpha}(P)$, we see that the outlyingness functions are equivalently bounded, i.e.

[TABLE]

Consequently the outlyingness function $O(f,P)$ is also bounded

[TABLE]

hence $D_{a}(f,P)\geq\alpha$ and so $[q]$ is in $D_{a,\alpha}(P)$. Thus the level sets induced on the amplitudes are convex. A similar proof shows that phase level sets are convex as well because the phase distance is convex. ∎

Proof of Proposition 3.4 (Uniform Consistency)

Proof.

Define the $\epsilon$-bracket $[l,u]$ as the set of all functions $f\in F_{M}$ such that $l<f<u$ and $|O_{a}(u,P)-O_{a}(l,P)|\leq\epsilon$. Because $F_{M}$ is a set of differentiable functions, only a finite number of $\epsilon$-brackets are needed to cover $F_{M}$. Therefore for any $f\in F_{M}$ there exists bracket $[u_{i},l_{i}]$ such that

[TABLE]

Consequently,

[TABLE]

By the strong law for sample quantiles $O_{n}(f,\mathbbm{P})\xrightarrow{a.s.}O_{a}(f,\mathbbm{P})$ for any fixed $f\in F_{M}$ so the right hand side goes to $\epsilon$ almost surely. Therefore if we take an $\epsilon$ sequence converging to 0 we get $\sup_{f}\{O_{a,n}(f,P)-O_{a}(f,P)\}\xrightarrow{a.s.}0$. Since $|O_{a,n}(f,P)-O_{a}(f,P)|$ forms an upper bound on $|D_{a,n}(f,P_{n})-D_{a}(f,P)|$ we get the result

[TABLE]

Similarly for the phase depth $D_{p}(\cdot,P)$,

[TABLE]

∎

A.3 Contamination by a Single Anomaly

We first considered the case when a single outlier is present in the data. We compared each method on their ability to rank the outlier as the most outlying function. For each of the six outlier models, we sampled 99 functions from the main model and 1 function from the contamination model, with compositional noise and magnitude outliers added to both. Each trajectory was sampled on the same equidistant 30 point grid over $[0,1]$. 1000 simulations were used for each outlier model to estimate the average ranks in Table 1.

The results are sorted so that the methods with the lowest average ranking (ED-A) across all models are listed first and those with the highest are last. On each of the seven outlier models, except for model 7, ED-A had average rankings very near 1.000. This means that in most scenarios the amplitude depths correctly rank the outlier as the most outlying function. Other methods such as TVD and ST-T1 succeed on models 1, 5, and 6 where they almost always correctly rank the outlier but fail on models 2 and 4 where their average ranks are far from 1.000.

The results on models 2 and 4 demonstrate an important feature of the elastic depths, which is not shared by other methods. Because the elastic depths are based on proper distance metrics in the amplitude and phase spaces, they can detect trajectories that lack amplitude and phase characteristics present in the rest of the data. Having either too low of an amplitude (phase) or too high of an amplitude (phase) is both indicative of shape outlyingness. Prior methods have primarily focused on trajectories with too high of an amplitude and thus struggle to adequately handle models 2 and 4.

A.4 Full Simulations

We now describe the full simulation setup used to compare the Elastic Depths. Total Variation Depth (TVD), the magnitude-shape plot (MS), the directional outlyingness measure (DIR), and the Geometric boxplots (GEOM) were implemented as described in Section 5.1. The three sequential transformations (ST-T1, ST-T2, and ST-D1) were implemented by hand and outliers were identified using functional boxplot function fbplot in the fda package with functional directional quantiles as the underlying depth measure. The functional outlier map (FOM) was created using the fom function from the mrfDepth package using default values and the functional directional quantiles as the underlying depth measure. The Robust Functional Tangential Angle Psuedo-depth (rFUNTA) was computed using the author’s bootstrapping procedure with $probs=0.01$. The order extended Integrated and Infimal depths (FDJ and IDJ respectively) were computed using the depthf.fd1 function with Half-Space (Tukey) depth from the ddalpha R package. Outliers were found using the authors described procedure with $order=2$ and a 1.5 IQR multiplier.

A.5 Parameter sensitivity analysis

The boxplot methodology proposed in Section 4.1 requires pre-specifying a multiplier $K$ on the IQR to determine the outlier cutoff value. In the traditional univariate case, this multiplier is set $1.5$ because, for Gaussian data, 1.5 times the IQR yields in interval covering approximately 99.3% of the data. Theoretical justification of a similar multiplier is much more difficult in the functional data case, so in this section we numerically investigate the robustness of our outlier detection method to different multipliers.

We first looked at the true positive rate of the elastic depths as a function of $k$. We have seen that any $k$ value less than around 2.25 (rightmost vertical black line) will have a nearly 100% true positive rate. This is because the outliers will generally receive much lower depth values than the inliers, thus the only scenario in which the outliers are misclassified as inliers is when $k$ is large so that the whisker surpasses the outliers.

We also looked at how the coverage of our boxplots changed as a function of $k$ (Figure 9). We ran the simulations again for the 6 amplitude models in the paper, except this time we allowed the boxplot multiplier $k$ to vary from 1 to 3. We measured the whisker’s coverage using the True Negative Rate (TNR), which is the percentage of inliers classified as inliers.

Overall, we found the coverage (TNR) of the whisker to be fairly insensitive to the value of $k$. Across each model, the TNR steadily rises from about TNR = 0.95 (slight under coverage) at $k=1$ until it asymptotes to TNR = 1 as $k$ approaches 3.

Models 1, 2, and 6 required larger values of $k$, than models 3, 4, and 5 to achieve the desired coverage of 99.3%. However, we found this mostly had to do with the sampling rate of the observed trajectories. Essentially the trajectories in Models 1, 2, and 6 were undersampled relative to their frequency content (complexity) which caused the depth distribution to become left skewed and the boxplot to undercover the distribution.

When trajectories are sampled below their Nyquist rate, for example, twice the maximum frequency component present in the trajectory, they are considered under-sampled and aliasing occurs. This leads to poor alignment, which has the net effect of skewing the amplitude depth distribution due to underestimating many of the shape distances. In such a case, as long as the sampling rate is increased, then the amplitude depth distribution will become more symmetric and a lower value of $k$ can be used to achieve 99.3% coverage. We demonstrate this in Figure 11.

A.6 Additional phase simulations

In Section 5.2 we noted that our phase simulations did not optimally represent pure phase differences. In fact, the phase outliers were more outlying in amplitude space than in phase space. To study pure phase differences better, we compared our method ED-P against TVD on a simple sin waves vs cosine waves example.

We generated inliers from a GP with a $\sin(2\pi t),t\in[0,1]$ mean function and outliers from a GP with a $\cos(2\pi t)$ mean function. We used the exponential covariance function with $r=0.5$ for the GPs covariance functions. 100 replications were performed to generate panel B of Figure 12. TVD’s detection skill deteriorated significantly from the results in the manuscript, while ED-P did not.

This result demonstrates the TVD was relying on non-phase differences to help differentiate phase outliers. Because, if TVD was not using non-phase information, then the loss of non-phase information would not have impacted its $F_{1}$ scores.

A.7 Examples

To demonstrate the real-world performance of elastic depth we applied it to three different data sets, one for each manifold. The first data set is a collection of U.S. Treasury yield curves ($\mathbbm{R}$), the second is a sample from the MPEG-7 image data set ($\mathbbm{R}^{2}$), and the last is hurricane trajectories across the Atlantic Ocean ($\mathbb{S}^{2}$). We show that the elastic depths can be applied consistently across each of the three manifolds and analyzed using the same boxplot methodology.

A.8 U.S Daily Treasury Yield Curve Rates ($\mathbbm{R}$)

We first consider the daily U.S. Treasury yield curves from January 2017 to April 2019 (United States Department of the Treasury,, 2019). The daily yield curve is a plot of bond terms, or time to maturity, against the associated interest rate on a given day. The shape of the yield curve has long been taken as an indicator of economic activity. Ordinarily, the yield curve is monotonically increasing as a function of time to maturity. However, in the months preceding recessions, the yield curve often becomes “flattened” and then “inverted” meaning that bonds with shorter maturity dates start to command higher interest rates than bonds with longer maturity dates.

We collected yield curve data for each day between January 1st, 2017 to April 2019 with terms spanning from 1 month to 30 years. Though the yield curve is only supported on a finite set of points, we treat it as though it were a continuous trajectory supported on a compact subset of $\mathbbm{R}$. This assumption is not altogether unreasonable given the relatively smooth relationship between time to maturity and interest rate. We then applied the elastic depth based boxplots with no thresholding to perform outlier detection, see Figure 13.

There are two sets of outliers identified in Figure 13. The first is the yield curve on 09/05/17 where the 1 Month interest rate spiked over the 3 Month interest rate. Though the spike was enough to force that day’s yield curve to become outlying, the overall shape of the curve is still monotonic, so this type of outlier is unlikely to be of practical interest. The second set of outliers is the group of 26 yield curves, also highlighted in red, that corresponds to the end of the observation period. Though these curves are not complete inversions they are considered “flat” since the 3 Month and 30 Year interest rates are similar. Because these flat curves follow a period of regular yield curve behavior they may be taken as a sign of an oncoming yield curve inversion and potentially an economic recession.

A.9 MPEG-7 Shape Data ($\mathbbm{R}^{2}$)

Our next data example comes from the MPEG-7 shape data set (Manjunath et al.,, 2002). This data set consists of 1300 trajectories in $\mathbbm{R}^{2}$, corresponding to 65 shape classes with 20 observations each. Each trajectory is an outline of some object such as an apple, turtle, or a butterfly sampled from frames in a video. Objects may be rotated, distorted, or magnified with respect to other objects of the same class.

To illustrate our method, we selected 16 shape classes from the available 65. We used a depth threshold of 0.05 so that at most one outlier would be detected within each shape class. The identified amplitude outliers are displayed in red against the inlying shapes in blue on the left hand side of Figure 14. We can see that, depending on the class, elastic depth can identify anomalous trajectories with both too high and too low amplitudes with respect to the rest of the set. For instance, the outlying jellyfish in the second column of row four has its tentacles spreading out from all sides of its body, whereas the inlying jellyfish all have their tentacles on the right. On the other hand, the star in the fourth column of row two lacks the larger rounded features of the inlying set which is why it too was identified as outlying.

In some of the classes, such as the bone or the pentagon, it is not visually obvious why the identified outliers are outlying in shape. To understand why these trajectories may be shaped differently than the rest, we plotted the most outlying trajectory (red) against the most central trajectory (blue) for each class in the right hand side of Figure 14. This contrast makes clear that the bone outlier lacks a protrusion on the top and the pentagon outlier lacks spokes, features which would make them far in amplitude from other trajectories in their class. In other classes, such as the butterfly, there is no single feature that distinguishes the outliers so its outlyingness results from the culmination of numerous small shape differences.