Regularized Halfspace Depth for Functional Data

In the last several decades, data depth methods for multivariate data have been developed and are shown to be a powerful tool for both exploratory data analysis and robust statistics. Many different depth notions have been proposed, such as Tukey’s halfspace depth (Tukey, 1975), simplicial depth (Liu, 1990), and Mahalanobis depth (Liu, 1992). These multivariate depth functions have not only been well-investigated theoretically (Donoho and Gasko, 1992, Zuo and Serfling, 2000) but have also been applied for handling different statistical problems (Yeh and Singh, 1997, Rousseeuw et al., 1999, Li et al., 2012).

In particular, Tukey’s halfspace depth has been popular due to its many desirable properties (Zuo and Serfling, 2000) and the robustness of its depth median (Donoho and Gasko, 1992). However, its generalization to functional data inevitably encounters a degeneracy issue, namely, a naive extension of Tukey’s depth has zero depth values almost surely. Intuitively, this is due to the infinite-dimensionality of the functional space, where the set of projections is “too rich” and thus Tukey’s depth defined from the most extreme projection is too small. A numerical illustration of the degeneracy phenomenon in practice is shown in Figure 1, where the naive direct extension of Tukey’s depth to functional data (red curve) assigns zero depth to the majority of observations, preventing depth comparisons. The degeneracy was first discovered by Dutta et al. (2011) and further investigated by Kuelbs and Zinn (2013), Chakraborty and Chaudhuri (2014a). See Table 1 of Gijbels and Nagy (2017) for a summary of the degeneracy issues in this and other existing functional depth notions. Tukey’s halfspace depth has therefore been practically dismissed from consideration for functional data due to the failure of the naive generalization.

In this work, we propose a new notion of functional depth called the regularized halfspace depth (RHD). We resolve the degeneracy issue by restricting the set of projection directions in the definition of the halfspace depth to the set of directions with reproducing kernel Hilbert space (RKHS) norm less than a given positive number $\lambda\in(0,\infty)$. Indeed, the proposed RHD obtained in this fashion is non-degenerate over a general class of random functions. This allows for depth-based rankings if, for example, the underlying distribution of the random element is Gaussian. Figure 1 illustrates that the RHD is free of the degeneracy issue as shown by the green and yellow curves for small and large $\lambda$ values; details about this simulation setup can be found in Section S2.1 in the Supplemental Materials. Smaller $\lambda$ results in greater regularization, in which case the distribution of the depth values become more positive and dispersed, while the original Tukey’s halfspace depth, corresponding to the red curve with $\lambda=\infty$, exhibits a degenerate behavior, since almost all observations have zero depth.

A related but highly different projection-based Tukey’s depth for functional data has been considered by Cuesta-Albertos and Nieto-Reyes (2008), who proposed the random Tukey depth which defines the depth as the most severe extremity as seen along a fixed number of random projections. The number of projections should not diverge to infinity because otherwise the random Tukey depth will approach Tukey’s halfspace depth, which degenerates given functional observations. In contrast, our approach targets to utilize infinitely many projections and thus is more sensitive against different types of extremeness; the degeneracy issue is resolved through regularizing the projections. These approaches both demonstrate the necessity of applying restrictions when defining functional depth through the infimum over projection directions.

In addition to non-degeneracy, we establish other desirable depth properties for the RHD in analogy to the ones for multivariate and functional data postulated by Zuo and Serfling (2000), Nieto-Reyes and Battey (2016), and Gijbels and Nagy (2017), respectively. We also include topological properties of the depth level sets. A uniform consistency of the sample RHD for the population RHD is verified over a totally bounded subset. A version of the Glivenko–Cantelli theorem is derived for a class of infinite-dimensional halfspaces as a by-product. In all, our new depth satisfies all of the desirable properties considered by Gijbels and Nagy (2017), unlike most other existing functional depths.

A prominent feature of the RHD is its flexibility in ranking curves based on different shape features controlled by the regularization parameter $\lambda$: RHD with a larger $\lambda$ emphasizes higher-frequency shape differences in the samples, while a smaller $\lambda$ emphasizes overall magnitude. The definition of “shape” is data-driven and comes from the modes of variation in the observed data. This feature leads to our new definitions of extremeness for curves informed by data, which is utilized and demonstrated in other applications such as outlier detection method and will be introduced next.

In recent years, substantial development of outlier detection methods in functional data have been accomplished by using notions of functional depth. The functional boxplot proposed by Sun and Genton (2011) is an extension of the original boxplot (cf. Tukey, 1977), where functional central region and whiskers are constructed based on a functional depth method, in particular the modified band depth (MBD) (López-Pintado and Romo, 2009). Outlier detection methods for functions taking values in a multivariate space have also been considered, by e.g., Hubert et al. (2015).

Trajectories outlying in only shape but not at any given time point are challenging to detect; Nagy et al. (2017) focused on the detection of this type of outliers by applying data depths to the joint distribution of a set of two or three time points altogether. Dai et al. (2020) proposed to apply transformations sequentially to the functional data to reveal shape outliers that are not easily detected by other methods. Alternative outlier detection methods based on functional data depth include functional bagplots (Hyndman and Shang, 2010), the outliergram (Arribas-Gil and Romo, 2014), and directional outlyingness (Dai and Genton, 2018).

Unlike magnitude outliers that tend to have an extreme average value over the support, shape outliers are hard to define in the first place because there is vast possibility of outlyingness in shape, in terms of, e.g., jumps, peaks, smoothness, trends etc. A shape outlier can often be inlying pointwise but has a different overall pattern than the majority of the sample curves over the support. A major benefit of the proposed RHD is that it captures different modes of variation in the sample $\mathcal{X}_{n}\equiv\{X_{i}\}_{i=1}^{n}$ depending on the choice of the regularization parameter $\lambda$. We propose a data-adaptive definition of shape outlyingness and establish a novel shape outlier detection method based on the proposed RHD. The method uses a large number of randomly drawn directions to approximate the sample RHD and applies the univariate boxplot to the data functions projected onto those directions. Acknowledging the reviewers’ recommendations, we have also considered phase variation and phase outliers. Namely, the observed curves exhibit variability in the x-axis, for example, when each curve is shifted (or delayed) compared to the others. As provided in the numerical studies, the proposed outlier detection method well detects different types of extremeness/outlyingness in shape, phase, as well as in magnitude.

Our methods are shown to be practical and useful with different real data examples: medfly, world population growth, and sea surface temperature datasets. Due to the complexity of the first dataset, it has barely been studied (cf. Sguera and López-Pintado, 2021) in the literature of functional outlier detection, while the other two datasets have been analyzed a few times previously (cf. Sun and Genton, 2011, Nagy et al., 2017, Dai et al., 2020). Our data illustrations demonstrate that the proposed method is relevant for analyzing both complicated and rough trajectories as well as relatively simple smooth curves.

Next we highlight our main contributions while summarizing the organization of the paper. Section 2 defines and studies the proposed RHD, which is a novel, flexible and useful depth notion for functional data; it is shown in Section 2.1 that the proposed RHD is non-degenerate over a large class of random functions. We then establish desirable depth properties of the RHD in Sections 2.2-2.3 and provide a theoretically reasonable and practical algorithm to compute it in Section 2.4; some intermediate results such as a Glivenko-Cantelli theorem are new in the functional depth literature. The choice of relevant tuning parameters are discussed in Section 2.5. Section 3 is devoted to describing the proposed outlier detection method; the details are presented in Section 3.1 along with a rule for selecting the adjustment factor. Section 3.2 emphasizes the robustness of the overall procedure against outliers. The proposed methods are useful and outstanding for identifying functional shape outliers, which can be challenging for the existing depth notions, as illustrated in the numerical study in Section 4. Data applications are described in Section 5, indicating that our methods are applicable to some complex structured functional data. All proofs, additional simulations, and extra results from the real data analyses can be found in the supplement.

Let $X$ be a random function defined on a probability space $(\Omega,\mathcal{F},\mathsf{P})$ that takes values in the underlying Hilbert space $\mathbb{H}$ with inner product $\langle\cdot,\cdot\rangle$. We suppose that $\mathbb{H}$ is infinite-dimensional and separable. The original Tukey’s halfspace depth (Tukey, 1975) of $x\in\mathbb{H}$ with respect to the probability measure $P_{X}$ induced by $X$ is, for $x\in\mathbb{H}$,

However, the direction set $\mathbb{S}\equiv\{v\in\mathbb{H}:\|v\|=1\}$ in the definition of the Tukey’s halfspace depth (1) is too rich due to the infinite dimension of $\mathbb{H}$. This set $\mathbb{S}$ contains an infinite orthonormal set in $\mathbb{H}$, which results in $\mathbb{S}$ not being totally bounded (MacCluer, 2009, Chapter 4) and the infimum in (1) being too small. As a result, the naive definition in (1) may become degenerate in an infinite dimensional space (Dutta et al., 2011), as described next.

We overcome this degeneracy issue and define a new notion of halfspace depth by introducing a regularization step. Let $\mathbb{H}$ be an infinite-dimensional Hilbert space, and define by $x\otimes y:\mathbb{H}\to\mathbb{H}$, the tensor product of two elements $x,y\in\mathbb{H}$, as a bounded linear operator $(x\otimes y)(z)=\langle z,x\rangle y$ for $z\in\mathbb{H}$. Under a finite second moment assumption $\mathsf{E}[\|X\|^{2}]<\infty$, the covariance operator of $X$ is defined as $\Gamma\equiv\mathsf{E}[(X-\mathsf{E}[X])\otimes(X-\mathsf{E}[X])]$. Since the covariance operator $\Gamma$ is self-adjoint, non-negative definite, and compact, by spectral decomposition, it admits the eigen-decomposition $\Gamma=\sum_{j=1}^{\infty}\gamma_{j}(\phi_{j}\otimes\phi_{j})$, where $(\lambda_{j},\phi_{j})$ is the $j$-th eigenvalue–eigenfunction pair of $\Gamma$ with the properties that $\gamma_{1}\geq\gamma_{2}\geq\cdots\geq 0$, $\gamma_{j}\to 0$ as $j\to\infty$, and the eigenfunctions $\{\phi_{j}\}_{j=1}^{\infty}$ form a complete orthonormal system for the closure of the image of $\Gamma$ (Hsing and Eubank, 2015). We further assume that the eigenvalues are positive and strictly decreasing, i.e., $\gamma_{1}>\gamma_{2}>\cdots>0$, in order to avoid non-interesting cases, such as where random elements supported only on a finite-dimensional subspace of $\mathbb{H}$.

The proposed regularized halfspace depth is defined next.

We refer to the norm $\|v\|_{\mathbb{H}(\Gamma)}\equiv\|\Gamma^{-1/2}v\|$ as the reproducing kernel Hilbert space (RKHS) norm because if $\mathbb{H}=L^{2}([0,1])$ and $X$ has a continuous covariance function $G(s,t)=\mathsf{cov}(X(s),X(t))$, $s,t\in[0,1]$, then $\mathbb{H}(\Gamma)\equiv\{v\in L^{2}([0,1]):\|\Gamma^{-1/2}v\|<\infty\}$ is an RKHS with reproducing kernel $G$ (cf. Wahba, 1973) and an infinite-dimensional dense subspace of $\mathbb{H}$. The definitions $\|v\|_{\mathbb{H}(\Gamma)}$ and $\mathbb{H}(\Gamma)$ are generally valid given any Hilbert space $\mathbb{H}$.

However, the RKHS $\mathbb{H}(\Gamma)$ is still large since it contains all eigenfunctions, leading to extremely small halfspace probabilities and subsequent degeneracy. For this reason, we regularize the RKHS norm $\|v\|_{\mathbb{H}(\Gamma)}$ of projection direction $v$ using a pre-specified number $\lambda\in(0,\infty)$. The resulting set of projection directions used to define the RHD $D_{\lambda}$ is written as

The set $\mathcal{V}_{\lambda}$ is then reasonably small in the sense that it is totally bounded (cf. Proposition S4 in the supplement).

The regularization is the key difference of our method from other existing depth notions. For a given $\lambda\in(0,\infty)$, we restrict the set of unit-norm projection directions used in the halfspace probabilities to $\mathcal{V}_{\lambda}$ in (3), which requires the RKHS norm of the projections to be at most $\lambda$. This will resolve the degeneracy issue as shown below in Theorem 2. To the best of our knowledge, regularization has never been used in the depth literature in the context of either finite or infinite dimensional spaces.

Without regularization, degeneracy occurs for the original Tukey’s depth (1) in $\mathbb{H}$ as explained in Theorem 1. In contrast, our projection set $\mathcal{V}_{\lambda}$ from (3) is moderately sized, since it contains a collection of infinite linear combinations of the eigenfunctions with the constraint that the coefficients of non-leading eigenfunctions must be small. For instance, let $a_{j}=[(\gamma_{j}2^{-j})/(\sum_{k}\gamma_{k}2^{-k})]^{1/2}$ for $j=1,2,\dots$. If $\lambda\geq(\gamma_{1}/2)^{-1}$, then $\sum_{j}a_{j}\phi_{j}\in\mathcal{V}_{\lambda}$. However, $\phi_{j}\in\mathcal{V}_{\lambda}$ if and only if $\gamma_{j}^{-1/2}\leq\lambda$, so $\mathcal{V}_{\lambda}$ contains only finitely many eigenfunctions; see also Proposition S1 in the supplement. This enables/allows our proposed RHD to avoid the degeneracy issue. We will further show in Section 2.2 that RHD is a flexible data depth method with desirable theoretical properties.

Theorem 2 shows that the RHD is free of degeneracy. To state the theorem, let $\mu=\mathsf{E}[X]$ and $F_{v}$ denote the cumulative distribution of the standardized projection $\langle X-\mu,v\rangle/\|\Gamma^{1/2}v\|$. For each $\lambda,t\in(0,\infty)$, we consider the following condition on the tail probabilities of the standardized projections:

Condition $E(\lambda,t)$ for Theorem 2 ensures that the upper tail probabilities of the standardized projection are not too small across the direction set $\mathcal{V}_{\lambda}$. This tail probability condition is mild and general, which holds especially when the standardized projection distribution $F_{v}$ does not depend on $v\in\mathcal{V}_{\lambda}$. This encompasses a wide range of random elements, particularly those whose standardized projections follow a well-known distribution, those whose FPC scores are dependent, or those whose higher moments are infinite, as listed below.

In addition to resolving degeneracy, regularization parameter $\lambda$ also provides flexibility on emphasizing the shape and/or magnitude of the functional observations. Given a small $\lambda$, the RHD focuses more on the overall magnitude of the input functions since the projection set concentrates around the leading eigenfunctions, which tend to capture lower frequency variations. On the other hand, as $\lambda$ increases, the RHD increases its emphasis on variation along the higher-frequency directions given by the non-leading eigenfunctions. This property of the RHD will be utilized later when we develop a new method for shape outlier detection. The regularization parameter $\lambda\in(0,\infty)$ is a user-specific parameter that depends on the purposes of the analysis. The selection of $\lambda$ will be discussed in Section 2.5.

Surprisingly, the RHD is shown to be robust against outliers despite utilizing the covariance, which is generally assumed to exist in functional data literature but typically considered non-robust. Moreover, the RHD can be made more robust by restricting the directions in $\mathcal{V}\equiv\{v\in\mathbb{H}:\|v\|=1\}$ to a general subset of $\mathcal{V}$ that does not require the finite moments. One option is to substitute the RHKS norm $\|\Gamma^{-1/2}v\|$ with another RKHS norm with a different reproducing kernel. However, such general subsets may not reflect the variability in shape and magnitude of the functional observations. We show that by imposing a second moment condition, the RHD offers numerous practical benefits, such as detecting many types of shape outliers in complex functional data sets; see Section 4 for numerical performance. These advantages are usually not obtained when using a single functional depth such as modified band depth (López-Pintado and Romo, 2009), extremal depth (Narisetty and Nair, 2016), integrated depth (Nagy et al., 2017), and so on; nonetheless, using a combination of multiple depth functions can give similar performance but adds complexity (Nagy et al., 2017, Dai et al., 2020). Despite the covariance being non-robust, the proposed methods based on RHD exhibit resistance to many types of outliers and are shown to be useful in practice. This has been illustrated and discussed in Section 4.4 and Section 3.2, respectively.

Our proposed RHD can be constructed whenever the underlying function space is a separable Hilbert space, for example, $L^{2}\equiv L^{2}([0,1])=\left\{f:[0,1]\to\mathbb{R}|\int_{0}^{1}f(t)^{2}dt<\infty\right\}$ with inner product $\langle f_{1},f_{2}\rangle\equiv\int_{0}^{1}f_{1}(t)f_{2}(t)dt$ for $f_{1},f_{2}\in L^{2}$. More generally, the RHD is also applicable to functional data taking values in a multivariate space (cf. Claeskens et al., 2014, Hubert et al., 2015, López-Pintado et al., 2014) given an appropriately defined inner product. For solidity, in this work, we focus on univariate curves in $\mathbb{H}=L^{2}$.

For a depth function to be intuitively sound, it should provide a center-outward ordering of the data and satisfy some desirable properties. Liu (1990) discussed some properties satisfied by the proposed simplicial depth. Zuo and Serfling (2000) formally postulated a list of desirable properties that should be satisfied by a well-designed depth function and showed that the original Tukey’s halfspace depth satisfies all these properties. Nieto-Reyes and Battey (2016), Gijbels and Nagy (2017) among others extended such depth properties for functional data. We show some basic properties of the RHD in the following theorems which extend the postulated multivariate depth properties introduced by Zuo and Serfling (2000) to the infinite-dimensional Hilbert space $\mathbb{H}$. The proofs are deferred to Section S1.4 in the Supplementary Materials.

Theorem 3(a) extends the affine invariance property (cf. Zuo and Serfling, 2000) to an infinite dimensional space and establishes this property for the proposed RHD because any surjective isometry between two normed spaces is affine by the Mazui–Ulam theorem (Väisälä, 2003). Properties (b) and (c) in Theorem 3 respectively states that the RHD is maximized at the unique center $\mu$ and it is nonincreasing along any ray leaving from the deepest curve $\mu\in\mathbb{H}$. These properties are analogous to their Euclidean versions in Zuo and Serfling (2000), but the proof for (b) is more involved and uses properties of the projection set $\mathcal{V}_{\lambda}$. Theorem 3(d) is a weakened version of vanishing at infinity property where we need some restrictions such as the unnormalized direction set $\mathcal{U}_{\lambda}$ to limit the behavior of the depth function; see Section S1.4 in the Supplementary Materials for another version of the vanishing at infinity property with the finite-dimensional assumption and a counterexample.

The next theorem establishes the continuity of the depth function and topological properties of its level sets. For $\tau\in[0,\infty)$, we define the upper level set $G_{\tau,\lambda}$ of $D_{\lambda}$ at level $\tau$ as $G_{\tau,\lambda}\equiv\{x\in\mathbb{H}:D_{\lambda}(x)\geq\tau\}.$ The following theorem extends analogous properties for the finite-dimensional data case shown in Theorem 2.11 of Zuo and Serfling (2000).

The upper semi-continuity of the RHD $D_{\lambda}$ and the closedness (and hence, the completeness) of the level set $G_{\tau,\lambda}$ are closely related. The compactness follows from vanishing at infinity and the restriction onto a finite dimensional subspace of $\mathbb{H}$. Without this restriction, the compactness of the level sets may not be achieved as described in Section S1.4 of the Supplementary Materials since vanishing at infinity fails.

The population version of the RHD from (2) cannot be computed in practice. To define its sample version estimated from the observed functions, let $X_{1},\dots,X_{n}$ be independently and identically distributed (iid) copies of $X$ with empirical distribution $\hat{P}_{n}$ defined as $\hat{P}_{n}(B)\equiv n^{-1}\sum_{i=1}^{n}\mathbb{I}(X_{i}\in B)$ for a Borel set $B\in\mathcal{B}(\mathbb{H})$, where $\mathbb{I}$ denotes the indicator function.. The sample covariance operator $\hat{\Gamma}_{n}$ is defined as $\hat{\Gamma}_{n}\equiv n^{-1}\sum_{i=1}^{n}(X_{i}-\bar{X})\otimes(X_{i}-\bar{X})$ with sample mean $\bar{X}\equiv n^{-1}\sum_{i=1}^{n}X_{i}$. Write the spectral decomposition of the operator $\hat{\Gamma}_{n}$ as $\hat{\Gamma}_{n}=\sum_{j=1}^{n}\hat{\gamma}_{j}(\hat{\phi}_{j}\otimes\hat{\phi}_{j})$, where $(\hat{\lambda}_{j},\hat{\phi}_{j})$ is the $j$-th eigenpair of $\hat{\Gamma}_{n}$ for $j=1,\dots,n$, with $\hat{\gamma}_{1}\geq\cdots\geq\hat{\gamma}_{n}\geq 0$ (Hsing and Eubank, 2015).

To define the sample version of the RHD, the set of projection directions is obtained as

where $\hat{\Gamma}_{J}^{-1/2}\equiv\sum_{j=1}^{J}\hat{\gamma}_{j}^{-1/2}\hat{\pi}_{j}$ is a $J$-truncated inverse of the square root operator $\hat{\Gamma}_{n}^{1/2}=\sum_{j=1}^{n}\hat{\gamma}_{j}^{1/2}\hat{\pi}_{j}$. Here, the truncation parameter $J=J(n)$ may depend on the sample size $n$ with $J(n)\to\infty$ as $n\to\infty$. We define the sample RHD of $x\in\mathbb{H}$ (with respect to $\hat{P}_{n}$) as

To establish the uniform consistency of the sample RHD $\hat{D}_{\lambda,n}(x)$ for the population RHD $D_{\lambda}(x)$, we need the following continuity/anti-concentration assumption imposed on the distribution of the projection of $X$: for a regularization parameter $\lambda\in(0,\infty)$ and a totally bounded subset $\mathbb{B}\subseteq\mathbb{H}$,

This condition is satisfied by a large class of random functions $X$ where the projection scores $\langle X,\phi_{j}\rangle$, $j=1,2,\dots$ are independent and continuous; in particular, Gaussian random functions.

To verify the consistency theorem, in addition to Condition $C(\lambda,\mathbb{B})$, we impose some standard regularity conditions on $X$, denoted as Conditions (A1)-(A4) in Section S1.5.3 of the supplement. Similar or stronger conditions are commonly imposed in the functional principal component analysis literature for consistent estimation of the eigen-components (Cardot et al., 2007, Hall et al., 2007). These conditions are applied in the proofs involving perturbation theory for functional data (cf. Sections S1.5.3-S1.5.4 of the supplement). The consistency result is established in the following theorem.

The proof for this result is included in the supplement. A main component of the proof is to establish the Glivenko–Cantelli theorem for the empirical process $\{\hat{P}_{n}(H_{x,v}):v\in\mathcal{V}_{\lambda},x\in\mathbb{B}\}$ indexed by an infinite-dimensional class of halfspaces, where $H_{x,v}\equiv\{y\in\mathbb{H}:\langle y-x,v\rangle\geq 0\}$ denote the closed halfspace with normal vector $v\in\mathbb{H}$ containing $x\in\mathbb{H}$ on its boundary. Empirical process results for infinite-dimensional space are challenging to show except for the results on bracketing numbers of the classes of either smooth or monotone functions by van der Vaart and Wellner (1996, Section 2.7) and Glivenko–Cantelli and Donsker theorems over some finite dimensional subspace by Chakraborty and Chaudhuri (2014b). Developing empirical process theory for functional data may not be feasible without such restriction on the class of functions to reasonably small subclasses. The main reason is that, in the infinite-dimensional Hilbert space $\mathbb{H}$, sufficiently small balls (e.g., balls with radius smaller than $\sqrt{2}$) fail to cover the closed unit ball; this implies the closed unit ball in $\mathbb{H}$ is not only non-totally-bounded (e.g., MacCluer, 2009, Hsing and Eubank, 2015), but it may also not have a finite covering number. As a related discussion, we refer to Caramanis (2000), showing the closed unit ball in an infinite dimensional Banach space is not a Glivenko-Cantelli class. Here, we consider a totally bounded subset $\mathbb{B}$ for such restriction. The proof consists of showing the total boundedness of $\mathcal{V}_{\lambda}$ and asymptotic uniform equicontinuity of the stochastic process $\{\hat{P}_{n}(H_{x,v}):v\in\mathcal{V}_{\lambda},x\in\mathbb{B}\}$ (cf. Section S1.4.2 of the supplement), which is more involved than directly applying the covering or bracketing numbers (cf. van der Vaart and Wellner, 1996, Chapter 2). Please refer to the respective paragraphs following Theorems 5-6.

The sample RHD involves an infimum over the projection set $\hat{\mathcal{V}}_{\lambda,J}$ which makes the computation hard; even for the original Tukey’s depth in $\mathbb{R}^{d}$ the exact computation can be prohibitive for dimension $d\geq 4$ (Dyckerhoff and Mozharovskyi, 2016). To approximate the sample RHD, we adopt a random projection approach.

For $\hat{v}=\sum_{j=1}^{J}\hat{a}_{j}\hat{\phi}_{j}\in\mathrm{span}\{\hat{\phi}_{1},\dots,\hat{\phi}_{J}\}$ with $\hat{\bm{a}}=(\hat{a}_{1},\dots,\hat{a}_{J})^{\top}\in\mathbb{R}^{J}$, it holds that

where $\hat{\mathcal{R}}_{J}\equiv\mathrm{diag}(\hat{\gamma}_{1},\dots,\hat{\gamma}_{J})$. Then, $\hat{\mathcal{V}}_{\lambda,J}$ can be written as

where $\hat{\mathcal{A}}_{\lambda,J}\equiv\{\hat{\bm{a}}\in\mathbb{R}^{J}:\|\hat{\bm{a}}\|_{\mathbb{R}^{J}}=1,\|\hat{\mathcal{R}}_{J}^{-1/2}\hat{\bm{a}}\|_{\mathbb{R}^{J}}\leq\lambda\}$. Note that the set $\hat{\mathcal{A}}_{\lambda,J}$ is the intersection of the (surface of) unit sphere and a (solid) ellipsoidal region centered at the origin in $\mathbb{R}^{J}$. Then, the sample RHD is written as

where $\hat{\bm{X}}_{i}=[\langle X_{i},\hat{\phi}_{j}\rangle]_{1\leq j\leq J}$ and $\hat{\bm{x}}=[\langle x,\hat{\phi}_{j}\rangle]_{1\leq j\leq J}$ are $J$-dimensional vectors in $\mathbb{R}^{J}$. Although (6) is reminiscent of sample Tukey’s depth applied on the data $\hat{\bm{X}}_{i}$, they are remarkably different in their designs. The RHD targets infinite-dimensional data and allows for an increasing $J$. The depth values are insensitive to increase in $J$ when it is moderately large ($J\geq 6$ as demonstrated in our numerical studies) thanks to the regularization set which stabilizes the depth. On the other hand, Tukey’s depth is most often applied in a fixed low dimension; when the dimension increases, Tukey’s depth become less capable of distinguishing points because of the concentration of measure and the degeneracy issue.

To approximate the infimum, we draw $M$ random vectors $\hat{\bm{a}}_{m}=(\hat{a}_{m1},\dots,\hat{a}_{mJ})^{\top}\overset{\mathsf{iid}}{\sim}\nu$ for $m=1,\dots,M$ from some continuous distribution $\nu$ supported on $\hat{\mathcal{A}}_{\lambda,J}\subseteq\mathbb{R}^{J}$. The projection set is approximated by the finite set

Then, the sample depth $\hat{D}_{\lambda,n}(x)=\hat{D}_{\lambda,n}(x,\hat{P}_{n})$ in (6) is approximated by

Rejection sampling is applied to draw the projection directions $\hat{\bm{a}}_{m}$. To increase acceptance in the target region $\hat{\mathcal{A}}_{\lambda,J}$, we set the proposal distribution of projection directions $\hat{\bm{a}}_{m}$ to that of $\bm{z}/\|\bm{z}\|_{\mathbb{R}^{J}}$ where $\bm{z}\sim\mathsf{N}(\bm{0},\hat{\mathcal{R}}_{J})$, a non-isotropic distribution on the unit sphere.

The following theorem justifies the approximated sample RHD in (8) for a sufficiently large $M$, for which the proof is given in Section S1.6 in the supplement.

The key point of the proof is to show a type of asymptotic equicontinuity of a process $\{\hat{P}_{n}(H_{x,\hat{v}}):x\in\mathbb{B},\hat{v}\in\hat{\mathcal{V}}_{\lambda,J}\}$, where the index set also depends on the observed random functions $\mathcal{X}_{n}\equiv\{X_{i}\}_{i=1}^{n}$. Such processes with a changing index set following the sample size $n$ have rarely been studied. We derived this theorem by constructing a correspondence from $\hat{\mathcal{V}}_{\lambda,J}$ to $\mathcal{V}_{\lambda}$, relating this correspondence with the asymptotic equicontinuity of the original process $\{\hat{P}_{n}(H_{x,v}):x\in\mathbb{B},v\in\mathcal{V}_{\lambda}\}$, and using the total boundedness of $\mathcal{V}_{\lambda}$.

The population (2) and sample RHD (5) are monotonically decreasing as $\lambda$ increases. In order for the approximate RHD to enjoy the same property, namely $\lambda\mapsto\tilde{D}_{\lambda,n,M}(x)$ is decreasing in $\lambda$ when all other quantities are fixed, we employ the following procedure. Let $\Lambda\equiv\{\lambda_{1},\dots,\lambda_{Q}\}$ be a finite set of regularization parameters that we will use to evaluate the approximate RHD. We draw a single set of random directions $\bm{z}_{1},\dots,\bm{z}_{M}\sim\mathsf{N}(0,\hat{\mathcal{R}}_{J})$ and set $\hat{\bm{a}}_{m}=\bm{z}_{m}/\|\bm{z}_{m}\|_{\mathbb{R}^{J}}$ as explained above, and use the directions $\{\hat{\bm{a}}_{m}\}_{m=1}^{M}$ to compute all depth values $\{\tilde{D}_{\lambda,n,M}:\lambda\in\Lambda\}$. The monotonicity then follows from the definition of $\tilde{D}_{\lambda,n,M}$.

In practice, instead of specifying $\lambda$ for regularization, we use the $u$-quantile of the RKHS norms $\mathcal{N}_{M}\equiv\{\|\hat{\mathcal{R}}_{J}^{-1/2}\hat{\bm{a}}_{m}\|_{\mathbb{R}^{J}}\}_{m=1}^{M}$ as $\lambda$, where $u\in(0,1)$ is called the quantile level for regularization. The quantile level approach to determine $\lambda$ has several advantages; it depends on the data and avoids choosing values of $\lambda$ outsides of the range of the RKHS norms $\mathcal{N}_{M}$. In our numerical studies, we investigate the performance of RHD using different $u$-quantiles of $\mathcal{N}_{M}$ as regularization parameter $\lambda$.

Algorithm 1 summarizes the approximate computation.

The approximate sample RHD $\tilde{D}_{\lambda,n,M}$ (8) involves three tuning parameters: the regularization parameter $\lambda$ (or equivalently, the quantile level $u$), the truncation level $J$, and the number $M$ of random projection directions. The choice of the most influential parameter $\lambda\in(0,\infty)$ depends on the user’s analytic goals, as each $\lambda$ can provide different perspectives (e.g., magnitude versus shape) and the optimal $\lambda$ for one analysis (e.g., a rank test) may not be ideal for another (e.g., a classification method). Moreover, further insight about the data can be obtained by investigating a family of the RHDs $\{D_{\lambda}\}_{\lambda\in\Lambda}$ with multiple regularization parameters $\Lambda\subseteq(0,\infty)$. The less sensitive tuning parameter $J$ should be a sufficiently large integer, which may be determined by a common criterion such as the fraction of variance explained (FVE) (cf. Kokoszka and Reimherr, 2017, Section 12.2). Lastly, a reasonably large value for the final parameter $M$ ensures a sufficient approximation, as supported by Remark 3. We elaborate more details on the selection of these tuning parameters below.

The most important parameter, by far, is the regularization $\lambda$ since it determines how sensitive the depth is to different modes of variation in the sample, which is a useful feature for analysis. When $\lambda$ is small, projections are constrained to have large components only along the leading eigenfunctions, in which case the RHD emphasizes the overall magnitude rather than shape of the curves. In contrast, when $\lambda$ is large, the RHD considers projections along both leading and later eigenfunctions, so it increases its emphasis in shape relative to magnitude. This flexibility of the RHD is used and highlighted in our outlier detection method proposed in Section 3 and is demonstrated numerically in Section 4.3.

The choice of $\lambda$ may depend on the structure of a given dataset and the purpose of the study. If a dataset can be explained by a few modes of variation (say, mostly in terms of a vertical shift), a small regularization parameter can be chosen; otherwise a larger regularization parameter should be considered. If a class label response exists, the selection of $\lambda$ can be tied to depth-based classification performance (Ghosh and Chaudhuri, 2005, López-Pintado and Romo, 2006, Hubert et al., 2017). Parameter $\lambda$ could also be chosen by maximizing the power of RHD-based rank test (for general depth-based rank test, see Liu and Singh, 1993, López-Pintado and Romo, 2009) if training samples are available. In general, practitioners could inspect a range of $\lambda$ values and obtain different rankings and outlyingness from the depth values. This would help understand the dataset from different perspectives in terms of magnitude vs shape. For practical uses, we recommend as default to set $\lambda$ according to a quantile level $u=0.95$ of the RKHS norms $\mathcal{N}_{M}$, which is implemented in our R package. This choice tends to display extremeness in both magnitude and shape while avoiding degeneracy as shown in our numerical studies.

The selection of the truncation level $J$ is a secondary consideration. It should be chosen as a large enough value, e.g., based on the FVE, because $J$ is needed for approximating the population RKHS norm $\|\Gamma^{-1/2}v\|$ by the sample version $\|\hat{\Gamma}_{J}^{-1/2}v\|$. In Figure 2, we show the rankings of several types of outliers when using the approximate RHD based on different $\lambda$ and $J$ parameters to illustrate how sensitive the rankings are to the choice of these parameters. A simulated dataset containing non-Gaussian functional data is considered where one outlier is added at a time.

A magnitude and three shape outliers of different types were ranked, where the $x$-axis and different panels display varying $\lambda$ and $J$, respectively (see Section 4.1 for more details). We define the RHD-based (normalized) rank of a functional observation $X_{i}$ in the sample $\{X_{i}\}_{i=1}^{n}$ as $\tilde{R}_{i}/n$, where $\tilde{R}_{i}$ denotes the rank of $\tilde{D}_{\lambda,n,M}(X_{i})$ in $\{\tilde{D}_{\lambda,n,M}(X_{i})\}_{i=1}^{n}$, and therefore, higher ranks correspond to larger depth. Minimum ranks are assigned to curves with the same depths. We observe that the results are highly similar when $J$ is moderately large, for example, $J\geq 6$ in this case. In contrast, the rankings for the smallest and largest regularization $\lambda$ are markedly different; the magnitude outlier gets the lowest rank for all $\lambda$s but this is not the case for the shape outliers. The shape outliers are only detected with higher $\lambda$s because larger $\lambda$ emphasizes shape more strongly.

The least consequential parameter is $M$, which should be chosen to be as large as computation affords. In practice, for example, $M=1000$ is often enough to get a reasonable approximation (cf. Remark 3).

In sum, we recommend to choose a large $M$, a sufficiently large $J$ that may depend on the sample size $n$, and an appropriate $\lambda$ for the purposes of the analysis.

The detailed outlier detection procedure is described as follows. We start by defining the most extreme functions with the least (approximate) RHD values as outlier candidates, for which the index set is denoted as

Let $X_{i_{0}}$, $i_{0}\in\mathcal{I}_{\min,\lambda}$ be an outlier candidate in the sample $\mathcal{X}_{n}\equiv\{X_{i}\}_{i=1}^{n}$. To construct univariate boxplots, we keep track of only the most extreme projection directions for evaluating the halfspace probabilities at $x=X_{i_{0}}$ for outlier detection. Among directions in $\tilde{\mathcal{V}}_{\lambda,J,M}$ (7), suppose that there are $K=K_{i_{0}}$ directions, denoted as $\{v_{m_{k}}\}_{k=1}^{K}$, that minimize the halfspace probability when computing $\tilde{D}_{\lambda,n,M}(X_{i_{0}})$ in (8), that is,

For each direction $v_{m_{k}}$, we obtain projections $\{\langle X_{i},v_{m_{k}}\rangle\}_{i=1}^{n}$ and label outliers using a univariate boxplot with a given factor $f$ to construct the fence. The fence is constructed as $[Q_{1}-f\times IQR,Q_{3}+f\times IQR]$ where $Q_{1}$, $Q_{3}$, and IQR denoting the first quartile, third quartile, and the interquartile range of the projections, respectively. We label as outliers $\mathcal{O}_{i_{0},k}$ the functional observations with projection lying outside of the fence along this direction. The final set of all outliers is obtained as

The algorithm for outlier detection is given in Algorithm 2.