Nonparametric ridge estimation

Introduction

Multivariate data in many problems exhibit intrinsic lower dimensional structure. The existence of such structure is of great interest for dimension reduction, clustering and improved statistical inference, and the question of how to identify and characterize this structure is the focus of active research. A commonly used representation for low-dimensional structure is a smooth manifold. Unfortunately, estimating manifolds can be difficult even under mild assumptions. For instance, the rate of convergence for estimating a manifold with bounded curvature blurred by homogeneous Gaussian noise, is logarithmic [Genovese et al. (2012a)], meaning that an exponential amount of data are needed to attain a specified level of accuracy. In this paper, we offer a way to circumvent this problem. We define an object, which we call a hyper-ridge set that can be used to approximate the low-dimensional structure in a data set. We show that the hyper-ridge set captures the essential features of the underlying low-dimensional structure while being estimable from data at a polynomial rate.

Let $X_{1},\ldots,X_{n}$ be a sample from a probability density $p$ defined on an open subset of $D$-dimensional Euclidean space and let $\hat{p}$ be an estimate of the density. We will define hyper-ridge sets (called ridges for short) for both $p$ and $\hat{p}$, which we denote by $R$ and $\hat{R}$. We consider two cases that make different assumptions about $p$. In the hidden manifold case (see Figure 1), we assume that the density $p$ is derived by sampling from a $d<D$ dimensional manifold $M$ and adding $D$-dimensional noise. In the density ridge case, we look for ridges of a density without assuming any hidden manifold, simply as a way of finding structure in a point cloud, much like clustering. The goal in both cases is to estimate the hyper-ridge set. Although in the former case, we would ideally like to estimate $M$, this is not always feasible for reasonable sample sizes, so we use the ridge $R$ as a surrogate for $M$. We focus on estimating ridges from point cloud data; we do not consider image data in this paper.

A formal definition of a ridge is given in Section 2. Let $1\leq d<D$ be fixed. Loosely speaking, we define a $d$-dimensional hyper-ridge set of a density $p$ to be the points where the Hessian of $p$ has $D-d$ strongly negative eigenvalues and where the projection of the gradient on that subspace is zero. Put another way, the ridge is a local maximizer of the density when moving in the normal direction defined by the Hessian.

Yet another way to think about ridges is by analogy with modes. We can define a mode to be a point where the gradient is 0 and the second derivative is negative, that is, the eigenvalues of the Hessian are negative. The Hessian defines a $(D-d)$-dimensional normal space (corresponding to the $D-d$ smallest eigenvalues) and a $d$ dimensional tangent space. A ridge point has a projected gradient (the gradient in the direction of the normal) that is 0 and eigenvalues in the normal space that are negative. Modes are simply dimensional ridges.

Note that the density is not uniform over the ridge. Indeed, there can be modes (-dimensional ridges) within a ridge. What matters is that the function rises sharply as we approach the ridge (strongly negative eigenvalue).

One of the main points of this paper is that $R$ captures the essential features of $M$. If we can live with the slight bias in $R$, then it is better to estimate $R$ since $R$ can be estimated at a polynomial rate while $M$ can only be estimated at a logarithmic rate. Throughout this paper, we take the dimension of interest $d$ as fixed and given.

Many different and useful definitions of a “ridge” have been proposed; see the discussion of related work at the end of this section. We make no claim as to the uniqueness and optimality of ours. Our definition is motivated by four useful properties that we demonstrate in this paper: {longlist}[1.]

If $\hat{p}$ is close to $p$, then $\hat{R}$ is close to $R$ where $\hat{R}$ is the ridge of $\hat{p}$ and $R$ is the ridge of $p$.

If the data-generating distribution is concentrated near a manifold $M$, then the ridge $R$ approximates $M$ both geometrically and topologically.

$R$ can be estimated at a polynomial rate, even in cases where $M$ can be estimated at only a logarithmic rate.

The definition corresponds essentially with the algorithm derived by Ozertem and Erdogmus (2011). That is, our definition provides a mathematical formalization of their algorithm.

Our broad goal is to provide a theoretical framework for understanding the problem of estimating hyper-ridge sets. In particular, we show that the ridges of a kernel density estimator consistently estimate the ridges of the density, and we find and upper bound on the rate of convergence. The main results of this paper are (stated here informally):

Stability (Theorem 4). If two densities are sufficiently close together, their hyper-ridge sets are also close together.

Estimation (Theorem 5). There is an estimator $\hat{R}$ such that

where $\operatorname{\mathsf{Haus}}$ is the Hausdorff distance, defined in equation (9). Moreover, $\hat{R}$ is topologically similar to $R$ in the sense that small dilations of these sets are topologically similar.

Surrogate (Theorem 7). In the Hidden Manifold case with small noise variance $\sigma^{2}$ and assuming $M$ has no boundary, the hyper-ridge set of the density $p$ satisfies

and $R$ is topologically similar to $M$. Hence, when the noise $\sigma$ is small, the ridge is close to $M$. Note that we treat $M$ as fixed while $\sigma\to 0$. It then follows that

This leaves open the question of how to locate the ridges of the density estimator. Fortunately, this latter problem has recently been solved by Ozertem and Erdogmus (2011) who derived a practical algorithm called the subspace constrained mean shift (SCMS) algorithm for locating the ridges. Ozertem and Erdogmus (2011) derived their method assuming that the underlying density function is known (i.e., they did not discuss the effect of estimation error). We, instead, assume the density is estimated from a finite sample and adapt their algorithm accordingly by including a denoising step in which we discard points with low density. This paper provides a statistical justification for, and extension to, their algorithm. We introduce a modification of their algorithm called SuRF (Subspace Ridge Finder) that applies density estimation, followed by denoising, followed by SCMS.

Related work. Zero dimensional ridges are modes and in this case ridge finding reduces to mode estimation and SCMS reduces to the mean shift clustering algorithm [Fukunaga and Hostetler (1975), Cheng (1995), Li, Ray and Lindsay (2007), Chacón (2012)].

If the hidden structure is a manifold, then the process of finding the structure is known as manifold estimation or manifold learning. There is a large literature on manifold estimation and related techniques. Some useful references are Niyogi, Smale and Weinberger (2008) Caillerie et al. (2011), Genovese et al. (2009, 2012a, 2012b, 2012c), Tenenbaum, de Silva and Langford (2000), Roweis and Saul (2000) and references therein.

The notion of ridge finding spans many fields. Previous work on ridge finding in the statistics literature includes Cheng, Hall and Hartigan (2004), Hall, Peng and Rau (2001), Wegman and Luo (2002), Wegman, Carr and Luo (1993) and Hall, Qian and Titterington (1992). These papers focus on visualization and exploratory analysis. An issue that has been discussed extensively in the applied math and computer science literature is how to define a ridge. A detailed history and taxonomy is given in the text by Eberly (1996). Two important classes of ridges are watershed ridges, which are global in nature, and height ridges, which are locally defined. There is some debate about the virtues of various definitions. See, for example, Norgard and Bremer (2012), Peikert, Günther and Weinkauf (2012). Related definitions also appear in the fluid dynamics literature [Schindler et al. (2012)] and astronomy [Aragón-Calvo et al. (2010), Sousbie et al. (2008)]. There is also a literature on Reeb graphs [Ge et al. (2011)] and metric graphs [Aanjaneya et al. (2012), Lecci, Rinaldo and Wasserman (2013)]. Metric graph methods are ideal for representing intersecting filamentary structure but are much more sensitive to noise than the methods in this paper. It is not our intent in this paper to argue that one particular definition of ridge is optimal for all purposes. Rather, we use a particular definition which is well suited for studying the statistical estimation of ridges.

More generally, there is a vast literature on hunting for structure in point clouds and analyzing the shapes of densities. Without attempting to be exhaustive, some representative work includes Davenport et al. (2010), Klemelä (2009), Adams, Atanasov and Carlsson (2011), Chazal et al. (2011), Bendich, Wang and Mukherjee (2012).

Throughout the paper, we use symbols like $C,C_{0},C_{1},c,c_{0},c_{1},\ldots$ to denote generic positive constants whose value may be different in different expressions.

Model and ridges

In this section, we describe our assumptions about the data and give a formal definition of hyper-ridge sets, which we call ridges from now on. Further properties of ridges are stated and proved in Section 4.

The data generating process under model (4) is equivalent to the following steps: {longlist}[1.]

Draw $B$ from a $\operatorname{Bernoulli}(\eta)$.

If $B=0$, draw $X$ from a uniform distribution on ${\mathcal{K}}$.

If $B=1$, let $X=Z+\sigma\varepsilon$ where $Z\sim W$ and $\varepsilon$ is additional noise. Points $X_{i}$ drawn from $\operatorname{Unif}({\mathcal{K}})$ represent background clutter. Points $X_{i}$ drawn from $W\star\Phi_{\sigma}$ are noisy observations from $M$. When $M$ consists of a finite set of points, this can be thought of as a clustering model.

denote the eigenvalues of $H(x)$ and let $\Lambda(x)$ be the diagonal matrix whose diagonal elements are the eigenvalues. Write the spectral decomposition of $H(x)$ as $H(x)=U(x)\Lambda(x)U(x)^{T}$. Let $V(x)$ be the last $D-d$ columns of $U(x)$ (i.e., the columns corresponding to the $D-d$ smallest eigenvalues). If we write $U(x)=[V_{\diamond}(x)\dvtx V(x)]$ then we can write $H(x)=[V_{\diamond}(x)\dvtx V(x)]\Lambda(x)[V_{\diamond}(x)\dvtx\break V(x)]^{T}$. Let $L(x)\equiv L(H(x))=V(x)V(x)^{T}$ be the projector onto the linear space defined by the columns of $V(x)$. We call this the local normal space and the space spanned by $L^{\perp}(x)=I-L(x)=V_{\diamond}(x)V_{\diamond}(x)^{T}$ is the local tangent space. Define the projected gradient

If the vector field $G(x)$ is Lipschitz then by Theorem 3.39 of Irwin (1980), $G$ defines a global flow as follows. The flow is a family of functions $\phi(x,t)$ such that $\phi(x,0)=x$ and $\phi^{\prime}(x,0)=G(x)$ and $\phi(x,s+t)=\phi(\phi(x,t),s)$. The flow lines, or integral curves, partition the space (see Lemma 2) and at each $x$ where $G(x)$ is nonnull, there is a unique integral curve passing through $x$. Thus, there is one and only one flow line through each nonridge point. The intuition is that the flow passing through $x$ is a gradient ascent path moving toward higher values of $p$. Unlike the paths defined by the gradient $g$ which move toward modes, the paths defined by the projected gradient $G$ move toward ridges. The SCMS algorithm, which we describe later, can be thought of as approximating the flow with discrete, linear steps $x_{k+1}\leftarrow x_{k}+hG(x_{k})$. [A proof that the linear interpolation of these points approximates the flow in the case $d=0$ is given in Arias-Castro, Mason and Pelletier (2013).]

Definition: The ridge $R$ of dimension $d$ is given by $R=\{x\dvtx\|G(x)\|=0,\lambda_{d+1}(x)<0\}$.

Note that the ridge consists of the destinations of the integral curves: $y\in R$ if $\lim_{t\to\infty}\pi(t)=y$ for some $\pi$ satisfying (7).

Our definition is motivated by Ozertem and Erdogmus (2011) but is slightly different. They first define the $d$-critical points as those for which $\|G(x)\|=0$. They call a critical point regular if it is $d$-critical but not $(d-1)$-critical. Thus, a mode within a one-dimensional ridge is not regular. A regular point with $\lambda_{d+1}<0$ is called a principal point. According to our definition, the ridge lies between the critical set and the principal set. Thus, if a mode lies on a one-dimensional ridge, we include that point as part of the ridge.

2 Assumptions

We now record the main assumptions about the ridges that we will require for the results.

Assumption (A0) differentiability. For all $x$, $g(x)$, $H(x)$ and $H^{\prime}(x)$ exist.

Assumption (A1) eigengap. Let $B_{D}(x,\delta)$ denote a $D$-dimensional ball of radius $\delta$ centered at $x$ and let $R\oplus\delta=\bigcup_{x\in R}B_{D}(x,\delta)$. We assume that there exists $\beta>0$ and $\delta>0$ such that, for all $x\in R\oplus\delta$, $\lambda_{d+1}(x)<-\beta$ and $\lambda_{d}(x)-\lambda_{d+1}(x)>\beta$.

Assumption (A2) path smoothness. For each $x\in R\oplus\delta$,

Condition (A1) says that $p$ is sharply curved around the ridge in the $D-d$ dimensional space normal to the ridge. To give more intuition about the condition, consider the problem of estimating a mode in one dimension. At a mode $x$, we have that $p^{\prime}(x)=0$ and $p^{\prime\prime}(x)<0$. However, the mode cannot be uniformly consistently estimated by only requiring the second derivative to be negative since $p^{\prime\prime}(x)$ could be arbitrarily close to 0. Instead, one needs to assume that $p^{\prime\prime}(x)<-\beta$ for some positive constant $\beta$. Condition (A1) may be thought of as the analogous condition for a ridge. (A2) is a third derivative condition which implies that the paths cannot be too wiggly. (A2) also constrains the gradient from being too steep in the perpendicular direction. Note that these conditions are local: they hold in a size $\delta$ neighborhood around the ridge.

Technical background

Now we review some background. We recommend that the reader quickly skim this section and then refer back to it as needed.

where $\|\cdot\|$ is the Euclidean norm. Given two sets $A$ and $B$, the Hausdorff distance between $A$ and $B$ is

is called the $\varepsilon$-dilation of $A$. The dilation can be thought of as a smoothed version of $A$. For example, if there are any small holes in $A$, these will be filled in by forming the dilation $A\oplus\varepsilon$.

We use Hausdorff distance to measure the distance between sets for several reasons: it is the most commonly used distance between sets, it is a very strict distance and is analogous to the familiar $L_{\infty}$ distance between functions for sets.

2 Topological concepts

This subsection follows Chazal, Cohen-Steiner and Lieutier (2009) and Chazal and Lieutier (2005). The reach of a set $K$, denoted by $\operatorname{\mathsf{reach}}(K)$, is the largest $r>0$ such that each point in $K\oplus r$ has a unique projection onto $K$. A set with positive reach is, in a sense, a smooth set without self-intersections.

Now we describe a generalization of reach called $\mu$-reach. The key point is simply that the $\mu$-reach is weaker than reach. The full details can be found in the aforementioned references. Let $A$ be a compact set. Following Chazal and Lieutier (2005) define the gradient $\nabla_{A}(x)$ of $d_{A}(x)$ to be the usual gradient function whenever this is well defined. However, there may be points $x$ at which $d_{A}$ is not differentiable in the usual sense. In that case, define the gradient as follows. For $x\in A$ define $\nabla_{A}(x)=0$ for all $x\in A$. For $x\notin A$, let $\Gamma(x)=\{y\in A\dvtx\|x-y\|=d_{A}(x)\}$. Let $\Theta(x)$ be the center of the unique smallest closed ball containing $\Gamma(x)$. Define $\nabla_{A}(x)=\frac{x-\Theta(x)}{d_{A}(x)}$.

The critical points are the points at which $\nabla_{A}(x)=0$. The weak feature size $\operatorname{\mathsf{wfs}}(A)$ is the distance from $A$ to its closest critical point. For $0<\mu<1$, the $\mu$-reach $\operatorname{\mathsf{reach}}_{\mu}(A)$ is $\operatorname{\mathsf{reach}}_{\mu}(A)=\inf\{d\dvtx\chi(d)<\mu\}$ where $\chi(d)=\inf\{\|\nabla_{A}(x)\|\dvtx d_{A}(x)=d\}$. It can be shown that $\operatorname{\mathsf{reach}}_{\mu}$ is nonincreasing in $\mu$, that $\operatorname{\mathsf{wfs}}(A)=\lim_{\mu\to 0}\operatorname{\mathsf{reach}}_{\mu}(A)$ and that $\operatorname{\mathsf{reach}}(A)=\lim_{\mu\to 1}\operatorname{\mathsf{reach}}_{\mu}(A)$.

As a simple example, a circle $C$ with radius $r$ has $\operatorname{\mathsf{reach}}(R)=r$. However, if we bend the circle slightly to create a corner, the reach is 0 but, provided the kink is not too extreme, the $\mu$-reach is still positive. As another example, a straight line as infinite reach. Now suppose we add a corner as in Figure 4. This set has 0 reach but has positive $\mu$-reach.

Two maps $f\dvtx A\to B$ and $g\dvtx A\to B$ are homotopic if there exists a continuous map $H\dvtx\times A\to B$ such that $H(0,x)=f(x)$ and $H(1,x)=g(x)$. Two sets $A$ and $B$ are homotopy equivalent if there are continuous maps $f\dvtx A\to B$ and $g\dvtx B\to A$ such that the following is true: (i) $g\circ f$ is homotopic to the identity map on $A$ and (ii) $f\circ g$ is homotopic to the identity map on $B$. In this case we write $A\cong B$. Sometimes $A$ fails to be homotopic to $B$ but $A$ is homotopic to $B\oplus\delta$ for every sufficiently small $\delta>0$. This happens because $B\oplus\delta$ is slightly smoother than $B$. If $A\cong B\oplus\delta$ for all small $\delta>0$, we will say that $A$ and $B$ are nearly homotopic and we will write A\approx^{{\mbox{\sim}}}B.

The following result [Theorem 4.6 in Chazal, Cohen-Steiner and Lieutier (2009)] says that if a set $K$ is smooth and $\widetilde{K}$ is close to $K$, then a smoothed version of $\widetilde{K}$ is nearly homotopy equivalent to $K$.

Let $K$ and $\widetilde{K}$ be compact sets and let $\varepsilon=\operatorname{\mathsf{Haus}}(\widetilde{K},K)$. If

then (\widetilde{K}\oplus\alpha)\approx^{{\mbox{\sim}}}K.

3 Matrix theory

We make extensive use of matrix theory as can be found in Stewart and Sun (1990), Bhatia (1997), Horn and Johnson (2013) and Magnus and Neudecker (1988).

The $\operatorname{\mathsf{vec}}$ operator converts a matrix into a vector by stacking the columns. Thus, if $A$ is $m\times n$ then $\operatorname{\mathsf{vec}}(A)$ is a vector of length $mn$. Conversely, given a vector $a$ of length $mn$, let $[[a]]$ denote the $m\times n$ matrix obtained by stacking $a$ columnwise into matrix form. We can think of $[[a]]$ as the “anti-vec” operator.

If $A$ is $m\times n$ and $B$ is $p\times q$ then the Kronecker $A\otimes B$ is the $mp\times nq$ matrix

If $A$ and $B$ have the same dimensions, then the Hadamard product $C=A\circ B$ is defined by $C_{jk}=A_{jk}B_{jk}$.

Also, if $A(x)=f(x)I$ then $A^{\prime}(x)=\operatorname{\mathsf{vec}}(I)\otimes(\nabla f(x))^{T}$ where $\nabla f$ denotes the gradient of $f$.

Let $H$ be a $D\times D$ square, symmetric matrix with eigenvalues $\lambda_{1}\geq\cdots\geq\lambda_{D}$. Let $\widetilde{H}$ be another square, symmetric matrix with eigenvalues $\widetilde{\lambda}_{1}\geq\cdots\geq\widetilde{\lambda}_{D}$. By Weyl’s theorem [Theorem 4.3.1 of Horn and Johnson (2013)], we have that

Properties of ridges

In this section, we examine some of the properties of ridges as they were defined in Section 2 and show that, under appropriate conditions, if two functions are close together then their ridges are close and are topologically similar.

It will be convenient to parameterize the gradient ascent paths by arclength. Thus, let $s\equiv s(t)$ be the arclength from $\pi(t)$ to $\pi(\infty)$:

Let $t\equiv t(s)$ denote the inverse of $s(t)$. Note that

which is a restatement of (7) in the arclength parameterization.

2 Differentials

We will need derivatives of $g$, $H$, and $L$. The derivative of $g$ is the Hessian $H$. Recall from (13) that $H^{\prime}(x)=\frac{d\operatorname{\mathsf{vec}}(H(x))}{dx^{T}}$. We also need derivatives along the curve $\gamma$. The derivative of a functions $f$ along $\gamma$ is

Thus, the derivative of the gradient $g$ along $\gamma$ is

We will also need the derivative of $H$ in the direction of a vector $z$ which we will denote by

where $H=H(\gamma(s))$ and $E=(d/ds)H(\gamma(s))=H^{\prime}(x;z)$ with $z=\gamma^{\prime}(s)$.

3 Uniqueness of the γ𝛾\gamma paths

Conditions (A0)–(A2) imply that, for each $x\in(R\oplus\delta)-R$, there is a unique path $\gamma$ passing through $x$.

We will show that the vector field $G(x)$ is Lipschitz over $R\oplus\delta$. The result then follows from Theorem 3.39 of Irwin (1980). Recall that $G=Lg$ and $g$ is differentiable. It suffices to show that $L$ is differentiable over $R\oplus\delta$. Now $L(x)=L(H(x))$. It may be shown that, as a function of $H$, $L$ is Frechet differentiable. And $H$ is differentiable by assumption. By the chain rule, $L$ is differentiable as a function of $x$. Indeed, $dL/dx$ is the $D^{2}\times D$ matrix whose $j$th column is $\operatorname{\mathsf{vec}}(L^{\dagger}E_{j})$ where $E_{j}=[[H^{\prime}e_{j}]]$, $L^{\dagger}$ denotes the Frechet derivative, and $e_{j}$ is the vector which is 1 in the $j$th coordinate and zero otherwise.

4 Quadratic behavior

Conditions (A1) and (A2) imply that the function $p$ has quadratic-like behavior near the ridges. This property is needed for establishing the convergence of ridge estimators. In this section, we formalize this notion of quadratic behavior. Give a path $\gamma$, define the function

Thus, $\xi$ is simply the drop in the function $p$ along the curve $\gamma$ as we move away from the ridge. We write $\xi_{x}(s)$ if we want to emphasize that $\xi$ corresponds to the path $\gamma_{x}$ passing through the point $x$. Since $\xi\dvtx[0,\infty)\to[0,\infty)$, we define its derivatives in the usual way, that is, $\xi^{\prime}(s)=d\xi(s)/ds$.

Suppose that (A0)–(A2) hold. For all $x\in R\oplus\delta$, the following are true: {longlist}[1.]

$\xi^{\prime}(s)=\|G(\gamma(s))\|$ and $\xi^{\prime}(0)=0$.

$\xi(s)\geq\frac{\beta}{4}\|\gamma(0)-\gamma(s)\|^{2}$.

1. The first condition $\xi(0)=0$ is immediate from the definition.

Since the projected gradient is 0 at the ridge, we have that $\xi^{\prime}(0)=0$.

3. Note that $(\xi^{\prime}(s))^{2}=\|G_{s}\|^{2}=G_{s}^{T}G_{s}=g_{s}^{T}L_{s}g_{s}\equiv a(s)$. Differentiating both sides of this equation, we have that $2\xi^{\prime}(s)\xi^{\prime\prime}(s)=a^{\prime}(s)$, and hence

Since $L_{s}L_{s}=L_{s}$ we have that $\accentset{\mbox{{\large.}}}{L}_{s}=L_{s}\accentset{\mbox{{\large.}}}{L_{s}}+\accentset{\mbox{{\large.}}}{L}_{s}L_{s}$, and hence

Recall that $\accentset{\mbox{{\large.}}}{g}_{s}=-\frac{H_{s}G_{s}}{\|G_{s}\|}$. Thus,

4. The first term in $\xi^{\prime\prime}(s)$ is $-\frac{G_{s}^{T}H_{s}G_{s}}{\|G_{s}\|^{2}}$. Since $G$ is in the column space of $V$, $G_{s}^{T}H_{s}G_{s}=G_{s}^{T}(V_{s}\Lambda_{s}V_{s}^{T})G_{s}$ where $\Lambda_{s}=\operatorname{diag}(\lambda_{d+1}(\gamma(s)),\ldots,\lambda_{D}(\gamma(s)))$. Hence, from (A1),

Now we bound the second term $\frac{g_{s}^{T}\accentset{\mbox{{\large.}}}{L}_{s}G_{s}}{\|G_{s}\|}$. Since $L_{s}+L^{\perp}_{s}=I$ and $L_{s}G_{s}=G_{s}$, we have $g_{s}^{T}\accentset{\mbox{{\large.}}}{L}_{s}G_{s}=g_{s}^{T}L_{s}\accentset{\mbox{{\large.}}}{L}_{s}G_{s}+g_{s}^{T}L_{s}^{\perp}\accentset{\mbox{{\large.}}}{L}_{s}G_{s}=g_{s}^{T}L_{s}\accentset{\mbox{{\large.}}}{L}_{s}L_{s}G_{s}+g_{s}^{T}L_{s}^{\perp}\accentset{\mbox{{\large.}}}{L}_{s}L_{s}G_{s}$. Now $|g_{s}^{T}L_{s}\accentset{\mbox{{\large.}}}{L}_{s}L_{s}G_{s}|=0$. To see this, note that $L_{s}L_{s}=L_{s}$ implies $L_{s}\accentset{\mbox{{\large.}}}{L}_{s}+\accentset{\mbox{{\large.}}}{L}_{s}L_{s}=\accentset{\mbox{{\large.}}}{L}_{s}$ implies $L_{s}\accentset{\mbox{{\large.}}}{L}_{s}L_{s}+\accentset{\mbox{{\large.}}}{L}_{s}L_{s}=\accentset{\mbox{{\large.}}}{L}_{s}L_{s}$ implies $L_{s}\accentset{\mbox{{\large.}}}{L}_{s}L_{s}=0$. To bound $g_{s}^{T}L_{s}^{\perp}\accentset{\mbox{{\large.}}}{L}_{s}L_{s}G_{s}$ we proceed as follows. Let $E=(d/ds)H(\pi(\gamma(s)))=H^{\prime}(x;z)$ with $z=\gamma^{\prime}(s)$. Then, from Davis–Kahan,

5 Stability of ridges

Suppose that (A0)–(A2) hold for $p$ and that (A0) holds for $\widetilde{p}$. Let $\psi=\max\{\varepsilon,\varepsilon^{\prime},\varepsilon^{\prime\prime}\}$ and let $\Psi=\max\{\varepsilon,\varepsilon^{\prime},\varepsilon^{\prime\prime},\varepsilon^{\prime\prime\prime}\}$. When $\Psi$ is sufficiently small:

(1) Conditions (A1) and (A2) hold for $\widetilde{p}$.

(2) We have: $\operatorname{\mathsf{Haus}}(R,\widetilde{R})\leq\frac{2C\psi}{\beta}$.

(3) If $\operatorname{\mathsf{reach}}_{\mu}(R)>0$ for some $\mu>0$, then \widetilde{R}\oplus\frac{4\psi}{\mu^{2}}\approx^{{\mbox{\sim}}}R.

It follows that, $\|L-\widetilde{L}\|\leq C\varepsilon^{\prime\prime}$ and $\sup_{x}\|G(x)-\widetilde{G}(x)\|\leq C\psi$.

Now let $\widetilde{x}\in\widetilde{R}$. Thus, $\|\widetilde{G}(\widetilde{x})\|=0$, and hence $\|G(\widetilde{x})\|\leq C\psi$. Let $\gamma$ be the path through $\widetilde{x}$ so that $\gamma(s)=\widetilde{x}$ for some $s$. Let $r=\gamma(0)\in R$. From part 2 of Lemma 3, note that $\xi^{\prime}(s)=\|G(\widetilde{x})\|$. We have

for some $u$ between and $s$. Since $\xi^{\prime}(0)=0$, from part 4 of Lemma 3, $C\psi\geq s\xi^{\prime\prime}(u)\geq\frac{s\beta}{2}$ and so $\|r-\widetilde{x}\|\leq s\leq\frac{2C\psi}{\beta}$. Thus, $d(\widetilde{x},R)\leq\|r-\widetilde{x}\|\leq 2C\psi/\beta$.

Now let $x\in R$. The same argument shows that $d(x,\widetilde{R})\leq 2C\psi/\beta$ since (A1) and (A2) hold for $\widetilde{p}$.

(3) Choose any fixed $\kappa>0$ such that $\kappa<\frac{\mu^{2}}{5\mu^{2}+12}$. When $\Psi$ is sufficiently small, $\Psi\leq\kappa\operatorname{\mathsf{reach}}_{\mu}(K)$. Then \widetilde{R}\oplus\frac{4\psi}{\mu^{2}}\approx^{{\mbox{\sim}}}R from Theorem 1.

Ridges of density estimators

Now we consider estimating the ridges in the density ridge case (no hidden manifold). Let $X_{1},\ldots,X_{n}\sim P$ where $P$ has density $p$ and let

We assume that all derivatives of $p$ up to and including fifth degree are bounded and continuous. We also assume the conditions on the kernel in Gine and Guillou (2002) which are satisfied by all the usual kernels. Results on $\|p(x)-\hat{p}_{h}(x)\|_{\infty}$ are given, for example, in Prakasa Rao (1983), Giné and Guillou (2002) and Yukich (1985). The results in those references imply that

For the derivatives, rates are proved in the sense of mean squared error by Chacón, Duong and Wand (2011). They can be proved in the $L_{\infty}$ norm using the same techniques as in Prakasa Rao (1983), Giné and Guillou (2002) and Yukich (1985). The rates are:

[See Arias-Castro, Mason and Pelletier (2013), e.g.] Let $\psi_{n}=(\frac{\log n}{n})^{{2}/{(D+8)}}$. Choosing $h\asymp\sqrt{\psi_{n}}$ we get that $\varepsilon\asymp\varepsilon^{\prime}\asymp\varepsilon^{\prime\prime}\asymp O_{P}(\psi_{n})\mbox{ and }\varepsilon^{\prime\prime\prime}=o_{P}(1)$. From Theorem 4 and the rates above we have the following.

Let $\hat{R}^{*}=\hat{R}\cap(R\oplus\delta)$. Under the assumptions above and assuming that (A1) and (A2) hold, we have, with $h\asymp\sqrt{\psi_{n}}$ that

If $\operatorname{\mathsf{reach}}_{\mu}(R)>0$ then \hat{R}^{*}\oplus O(\psi_{n})\approx^{{\mbox{\sim}}}R.

Let $h>0$ be fixed and let $\widetilde{\psi}_{n}=\sqrt{\log n/n}$. Let $\hat{R}^{*}=\hat{R}\cap(R\oplus\delta)$. Under the assumptions above and assuming that (A1) and (A2) hold for $R_{h}$ we have, that

If $\operatorname{\mathsf{reach}}_{\mu}(R_{h})>0$ then \hat{R}^{*}\oplus O(\widetilde{\psi}_{n})\approx^{{\mbox{\sim}}}R.

Ridges as surrogates for hidden manifolds

We want to show that the ridge of $p_{\sigma}$ is a surrogate for $M$. Specifically, we show that, as $\sigma$ gets small, there is a subset $R_{*}\subset R$ in a neighborhood of $M$ such that $\operatorname{\mathsf{Haus}}(M,R_{*})=O(\sigma^{2}\log(1/\sigma))$ and such that R_{*}\approx^{{\mbox{\sim}}}M. We assume that $\eta=1$ in what follows; the extension to $0<\eta<1$ is straightforward. We also assume that $M$ is a compact $d$-manifold with positive reach $\kappa$. We need to assume that $M$ has positive reach rather than just positive $\mu$-reach. The reason is that, when $M$ has positive reach, the measure $W$ induces a smooth distribution on the tangent space $T_{x}M$ for each $x\in M$. We need this property in our proofs but this property is lost if $M$ only has positive $\mu$-reach for some $\mu<1$ due to the presence of unsmooth features such as corners.

where $\phi_{\sigma}(u)=(2\pi)^{-D/2}\sigma^{-D}$ $\exp(-\frac{\|u\|^{2}}{2\sigma^{2}})$. Thus, $p_{\sigma}$ is a mixture of Gaussians. However, it is a rather unusual mixture; it is a singular mixture of Gaussians since the mixing distribution $W$ is supported on a lower dimensional manifold.

Let $T_{x}M$ be the tangent space Recall that the tangent space at a point $x$ is the linear space spanned by the derivative vectors of smooth curves on the manifold through that point. to $M$ at $x$ and let $T_{x}^{\perp}M$ be the normal space to $M$ at $x$. Define the fiber at $x\in M$ by $F_{x}=T_{x}^{\perp}M\cap B_{D}(x,r)$. A consequence of the fact that the reach $\kappa$ is positive and $M$ has no boundary is that, for any $0<r<\kappa$, $M\oplus r$ can be written as a disjoint union

Let $r_{\sigma}>0$ satisfy the following conditions:

Specifically, take $r_{\sigma}=\alpha\sigma$ for some $0<\alpha<1$. Fix any $A\geq 2$ and define

Suppose that $\kappa=\operatorname{\mathsf{reach}}(M)>0$. Let $R_{\sigma}$ be the ridge set of $p_{\sigma}$. Let $M_{\sigma}=M\oplus r_{\sigma}$ and $R^{*}_{\sigma}=R_{\sigma}\cap M_{\sigma}$. For all small $\sigma>0$: {longlist}[1.]

$R^{*}_{\sigma}$ satisfies (A1) and (A2) with $\beta=c\sigma^{-(D-d+2)}$ form some $c>0$.

$\operatorname{\mathsf{Haus}}(M,R^{*}_{\sigma})=O(K_{\sigma}^{2})$.

R^{*}_{\sigma}\oplus CK_{\sigma}^{2}\approx^{{\mbox{\sim}}}M. If $R_{\sigma}$ is instead taken to be the ridge set of $\log p_{\sigma}$ then the same results are true with $\beta=c\sigma^{-2}$ and $M_{\sigma}=M\oplus\kappa$.

Without the assumption that $M$ has no boundary, there would be boundary effects of order $K_{\sigma}$. That is, the Hausdorff distance behaves like $O(K_{\sigma})$ for points near the boundary and like $O(K_{\sigma}^{2})$ for points not near the boundary.

The theorem shows that in a neighborhood of the manifold, there is a well-defined ridge, that the ridge is close to the manifold and is nearly homotopic to the manifold. It is interesting to compare the above result to recent work on finite mixtures of Gaussians [Carreira-Perpinan and Williams (2003), Edelsbrunner, Fasy and Rote (2012)]. In those papers, it is shown that there can be fewer or more modes than the number of Gaussian components in a finite mixture. However, for small $\sigma$, it is easy to see that for each component of the mixture, there is a nearby mode. Moreover, the density will be highly curved at those modes. Theorem 7 can be thought of as a version of the latter two facts for the case of manifold mixtures.

The theorem refers to the ridges defined by $p_{\sigma}$ and the ridges defined by $\log p_{\sigma}$. Although the location of the ridge sets is the same for both cases, the behavior of the function around the ridges is different. There are several reasons we might want to use $\log p$ rather than $p$. First, when $p$ is Gaussian, the ridges of $\log p$ correspond to the usual principal components. Second, the surrogate theorem holds in an $O(1)$ neighborhood of $M$ for the log-density whereas it only holds in an $O(\sigma)$ neighborhood of $M$ for the density.

To prove the theorem, we need a preliminary result. Let

Given a point $x$ let $\hat{x}$ be its projection onto $M$. In what follows, if $T$ is a matrix, then an expression of the form $T+O(r_{n})$ is to be interpreted to mean $T+B_{n}$ where $B_{n}$ is a matrix whose entries are of order $O(r_{n})$. Let

For all $x\in M_{\sigma}$, {longlist}[1.]

$p_{\sigma}(x)=\phi_{\perp}(x-\hat{x})(1+O(\widetilde{\sigma}))$.

Let $p_{\sigma,B}(x)=\int_{M\cap B}\phi_{\sigma}(x-z)\,dW(Z)$. Then $p_{\sigma,B}(x)=\phi_{\perp}(x-\hat{x})(1+O(\widetilde{\sigma}))$.

$g_{\sigma}(x)=-\frac{1}{\sigma^{2}}p_{\sigma}(x)((x-\hat{x})+O(K_{\sigma}^{2}))$ and $\|g_{\sigma}(x)\|=O(\sigma^{-(D-d-1)})$.

The projection matrix $L_{\sigma}$ satisfies

where $O_{\perp}(K^{2}_{\sigma})$ is a term of size $O(K_{\sigma}^{2})$ in $T_{x}^{\perp}$.

The proof is quite long and technical and so we relegate it to the Appendix.

Proof of Theorem 7 Let us begin with the ridge based on $p_{\sigma}$.

(1) Condition (A1) follows from parts 8 and 1 of Lemma 8 together with equation (49).

To verify (A2), we use parts 3 and 8 of Lemma 8: we get, for all small $\sigma$, that

(2) Suppose that $x\in R_{\sigma}^{*}$. Then $\|G_{\sigma}(x)\|=0$. Let $\hat{x}$ be the unique projection of $x$ onto $M$. From part 6 of Lemma 8,

Now let $\hat{x}\in M$. From the expression above, we see that $\|G_{\sigma}(\hat{x})\|=O(K_{\sigma}^{2})$. Let $\gamma$ be the path through $x$ and let $r$ be the destination of the path. Hence $\gamma(s)=x$ for some $s$ and $\gamma(0)=r$. Now we use Lemma 3. Then $\|G\|=\xi^{\prime}$ and

and so $\|x-\hat{r}\|=O(K_{\sigma}^{2})$. Hence, $\operatorname{\mathsf{Haus}}(R_{\sigma},M)=O(K_{\sigma}^{2})$.

(3) Homotopy. This follows from part (2) and Theorem 1.

Now consider the ridges of $\log p_{\sigma}(x)$. The proof is essentially the same as the proof above. The main difference is the Hessian as we now explain. Note that the Hessian $H_{\sigma}^{*}$ for $\log p_{\sigma}(x)$ is

From Lemma 8, parts 3 and 4, it follows that (after an appropriate rotation),

Notice in particular, that the dominant term of the smallest eigenvalue of $-\beta H_{\sigma}^{*}(x)$ is 1 whereas that the dominant term of the smallest eigenvalue of $-\beta H_{\sigma}(x)$ is 1 $d^{2}_{M}(x)/\sigma^{2}$ which is why we required $\|x-\hat{x}\|$ to be less than $\sigma$ in Theorem 7. Here, we only require that $\|x-\hat{x}\|\leq\kappa$.

We may now combine Theorems 4, 5, 6 and 7 to get the following.

Let $\hat{R}^{*}$ be defined as in Theorem 5. Then

Similarly, if $\hat{R}^{*}$ be defined as in Theorem 6 then

SuRFing the ridge

Here, we discuss Subspace Ridge Finding (SuRF) by using density estimation, followed by denoising and then followed by the subspace constrained mean shift (SCMS) algorithm due to Ozertem and Erdogmus (2011). We will not go into great details about the algorthm; we refer the reader to Ozertem and Erdogmus (2011).

Let us begin by reviewing the mean shift algorithm. The mean shift algorithm [Fukunaga and Hostetler (1975), Cheng (1995), Comaniciu and Meer (2002)] is a method for finding the modes of a density by approximating the steepest ascent paths. The algorithm starts with a mesh of points and then moves the points along the gradient ascent trajectories toward local maxima.

Given a sample $X_{1},\ldots,X_{n}$ from $p$, consider the kernel density estimator

where $K$ is a kernel and $h>0$ is a bandwidth. Let ${\mathcal{M}}=\{{v_{1},\ldots,v_{m}}\}$ be a collection of mesh points. These are often taken to be the same as the data but in general they need not be. Let $v_{j}(1)=v_{j}$ and for $t=1,2,3,\ldots$ we define the trajectory $v_{j}(1),v_{j}(2),\ldots,$ by

It can be shown that each trajectory $\{v_{j}(t)\dvtx t=1,2,3,\ldots,\}$ follows the gradient ascent path and converges to a mode of $\hat{p}_{h}$. Conversely, if the mesh ${\mathcal{M}}$ is rich enough, then for each mode of $\hat{p}_{h}$, some trajectory will converge to that mode.

The SCMS algorithm mimics the mean shift algorithm but it replaces the gradient with the projected gradient at each step. The algorithm can be applied to $\hat{p}$ or any monotone function of $\hat{p}$. As we explained earlier, there are some advantages to using $\log\hat{p}$. Figure 5 gives the algorithm for the log-density. This is the version we will use in our examples. Figure 6 gives the full SuRF algorithm.

The SCMS algorithm provides a numerical approximation to the paths $\gamma$ defined by the projected gradient. We illustrate the numerical algorithm in Section 8.

Implementation and examples

Here, we demonstrate ridge estimation in some two-dimensional examples. In each case, we will find the one-dimensional ridge set. Our purpose is to show proof of concept; there are many interesting implementation details that we will not address here. In each case, we use SuRF.

To implement the method requires that we choose a bandwidth $h$ for the kernel density estimator. There has been recent work on bandwidth selection for multivariate density estimators such as Chacón and Duong (2010, 2012) and Panaretos and Konis (2012). For the purposes of this paper, we simply use the Silverman rule [Scott (1992)].

Figures 7 through 10 show two examples of SuRF. In the first example, the manifold is a circle. Although the circle example may seem easy, we remind the reader that no existing statistical algorithms that we are aware of can, without prior assumptions, take a point cloud as input and find a circle, automatically.

The second example is a stylized “cosmic web” of intersecting line segments and with random background clutter. This is a difficult case that violates the assumptions; specifically the underlying object does not have positive reach. The starting points for the SCMS algorithm are a subset of the grid points at which a kernel density estimator is evaluated. We select those points for which the estimated density is above a threshold relative to the maximum value.

Figure 9 shows the estimator for four bandwidths. This shows an interesting phenomenon. When the bandwidth $h$ is large, the estimator is biased (as expected) but it is still homotopy equivalent to the true $M$. However, when $h$ gets too small, we see a phase transition where the estimator falls apart and degenerates into small pieces. This suggests it is safer to oversmooth and have a small amount of bias. The dangers of undersmoothing are greater than the dangers of oversmoothing.

The theory in Section 6 required the underlying structure to have positive reach which rules out intersections and corners. To see how the method fares when these assumptions are violated, see Figure 10. While the estimator is far from perfect, given the complexity of the example, the procedure does surprisingly well.

Conclusion

We presented an analysis of nonparametric ridge estimation. Our analysis had two main components: conditions that guarantee that the estimated ridge converges to the true ridge, and conditions to relate the ridge to an underlying hidden manifold.

We are currently investigating several questions. First, we are finding the minimax rate for this problem to establish whether or not our proposed method is optimal. Also, Klemelä (2005) has derived mode estimation procedures that adapt to the local regularity of the mode. It would be interesting to derive similar adaptive theory for ridges. Second, the hidden manifold case required that the manifold had positive reach. We are working on relaxing this condition to allow for corners and intersections (often known as stratified spaces). Third, we are developing an extension where ridges of each dimension $d=0,1,\ldots$ are found sequentially and removed one at a time. This leads to a decomposition of the point cloud into structures of increasing dimension. Finally, there are a number of methods for speeding up the mean shift algorithm. We are investigating how to adapt these speedups for SuRF.

As we mentioned in the Introduction, there is recent work on metric graph reconstruction which is a way of modeling intersecting filaments [Aanjaneya et al. (2012), Lecci, Rinaldo and Wasserman (2013)]. These algorithms have the advantage of being designed to handle intersecting ridges. However, it appears that they are very sensitive to noise. Currently, we are investigating the idea of first running SuRF and then applying metric graph reconstruction. Preliminary results suggest that this approach may get the best of both approaches.

Appendix

The purpose of this appendix is to prove Lemma 8. Recall that the gradient is $g_{\sigma}(x)=-\frac{1}{\sigma^{2}}\int_{M}(x-z)\phi_{\sigma}(x-z)\,dW(z)$ and the Hessian is

We can partition $M_{\sigma}$ into disjoint fibers. Choose an $x\in M_{\sigma}$ and let $\hat{x}$ be the unique projection of $x$ onto $M$. Let $B=B(\hat{x},K_{\sigma})$. For any bounded function $f(x,z)$,

Let $T=T_{\hat{x}}M$ denote the $d$-dimensional tangent space at $\hat{x}$ and let $T^{\perp}$ denote the $(D-d)$-dimensional normal space. For $z\in B\cap M$, let $\overline{z}$ be the projection of $z$ onto $T$. Then

where $u=(x-\hat{x})/d_{M}(x)\in T^{\perp}$ and $R=(\overline{z}-z)$. [Recall that $d_{M}$ is the distance function; see (8).] For small enough $\sigma$, there is a smooth map $h$ taking $z$ to $\overline{z}$ that is a bijection $B\cap M$ and so the distribution $W$ induces a distribution $\overline{W}$, that is, $\overline{W}(A)=W(h^{-1}(A))$. Let $\overline{w}$ denote the density of $\overline{W}$ with respect to Lebesgue measure $\mu_{d}$ on $T$. The density is bounded above and below and has two continuous derivatives.

For every $x\in R\oplus\sigma$, $\sup_{z\in B}\|z-\overline{z}\|\leq cK^{2}_{\sigma}$.

Recall that $r_{\sigma}=\alpha\sigma$ with $0<\alpha<1$. Define the following quantities:

First note that, for all $x\in R_{\sigma}$,

and so, $\phi_{\perp}(x-\hat{x})\asymp\sigma^{-(D-d)}$ as $\sigma\to 0$. Now,

Now $\|z-\overline{z}\|^{2}=O(K_{\sigma}^{4})$ and $|\langle x-\hat{x},\overline{z}-z\rangle|\leq\|x-\hat{x}\|\|\overline{z}-z\|=O(\sigma K^{2}_{\sigma})$ and so

Proof of Lemma 8 1. From (46), $p_{\sigma}(x)=\int_{M\cap B}\phi_{\sigma}(x-z)\,dW(z)+O(\sigma^{A})$. Now

2. $p_{\sigma,B}(x)$. This follows since in part 1 we showed that $p_{\sigma,B}(x)=p_{\sigma}(x)+O(\sigma^{A})$.

For some $u$ between $\hat{x}$ and $\overline{z,}$ we have

where $A=\{t=(\overline{z}-\hat{x})/\sigma\in h^{-1}(B)\}$. Finally,

It follow from part 1 that $\|g_{\sigma}(x)\|=O(\sigma^{-(D-d-1)})$.

4. To find the eigenvalues, we first approximate the Hessian. Without loss of generality, we can rotate the coordinates so that $T$ is spanned by $e_{1},\ldots,e_{d}$, $T^{\perp}$ is spanned by $e_{d+1},\ldots,e_{D}$ and $u=(0,\ldots,0,1)$. Now,

Let $Q=\int_{M\cap B}(x-z)(x-z)^{T}\phi_{\sigma}(x-z)\,dW(z)$. Then, from (47), we have $Q=Q_{1}+Q_{2}+Q_{3}+Q_{4}+Q_{5}+Q_{6}$ where

Next, with $t=(t_{1},\ldots,t_{d},0,\ldots,0)$,

A similar analysis on the remaining terms yields:

5. This follows from part 4 and the Davis–Kahan theorem.

6. From part 5, $L_{\sigma}(x)=L^{\dagger}+E$ where $L^{\dagger}=[{0_{d\times d}\enskip 0_{d,D-d}\atop 0_{D-d,d}\enskip I_{D-d}}]$ and $E=O(\widetilde{\sigma})$. Hence, $G_{\sigma}(x)=L_{\sigma}(x)g_{\sigma}(x)=(L^{\dagger}+E)g_{\sigma}(x)$ and the result follows from parts 3 and 4.

8. Now we turn to $\|H^{\prime}_{\sigma}\|$. Let $\Delta=(x-z)$. We claim that

To see this, note first that $H=\frac{1}{\sigma^{4}}Q-\frac{1}{\sigma^{2}}A$ where

Note that $Q=\int(x-z)(x-z)^{T}\Phi\,dW(z)$ where $\Phi=\phi_{\sigma}(\Delta)I_{D}$. So

Now $(d/dx)(\Delta\Delta^{T}\Phi)=(fg)^{\prime}$ where $f=\Delta\Delta^{T}$ and $g=\Phi$ and so

Each of these terms is of order $O(\sup_{x\in M}\|w^{\prime\prime}(x)\|/\sigma^{D-d+1})$. Consider the first term

where $u=(x-z)/\sigma$. As in the proof of part 1, we can restrict to $B\cap M$, do a change of measure to $\overline{W}$ and the term is dominated by

The other terms may be bounded similarly.

Acknowledgements

The authors thank the reviewers for many suggestions that improved the paper. In particular, we thank the Associate Editor who suggested a simplified proof of Lemma 3.

References