Stable Density Ridges: Consistency and Convergence of Subspace Constrained Mean Shift

Stable Density Ridges: Consistency and Convergence of Subspace Constrained Mean Shift Wanli Qiao

arXiv:2608.05112v1 [stat.ML] 5 Aug 2026

August 6, 2026 Abstract The Subspace Constrained Mean Shift (SCMS) algorithm is a popular nonparametric method for extracting density ridges, which serve as a low-dimensional representation of high-dimensional data. It is a widely held belief in the literature that SCMS trajectories converge to the classical density ridge, which we call the “static ridge”, defined via the density gradient and the eigenvalues and eigenvectors of the density’s Hessian. In this paper, we demonstrate that this assumption does not hold in general, as the static definition fails to account for the rotation of the trailing eigenspace along the continuous flow of the algorithm’s underlying vector field. To resolve this, we propose a paradigm shift by introducing the “stable ridge”, a novel geometric structure defined through the lens of dynamical systems and the Jacobian of the projected density gradient. We prove that this stable ridge is the true theoretical target of the SCMS algorithm. Building upon this foundation, we develop a generalized SCMS framework utilizing a constant step size, establishing its uniform R-linear convergence and topological surjectivity onto the stable ridge. We further derive the rates of convergence for estimating the stable ridge in terms of the Hausdorff distance. Finally, we expose that the original SCMS algorithm suffers from polynomial-time computational complexity, which is caused by implicitly coupling the step size to the smoothing bandwidth via the Mean Shift operator, and demonstrate how our generalized framework provides a statistically consistent and more efficient solution. Keywords: Subspace Constrained Mean Shift, Stable Density Ridges, Manifold Learning, Nonparametric Estimation, Computational Complexity, Dynamical Systems.

1

Introduction

The extraction of low-dimensional geometric structures from high-dimensional data is a fundamental problem in modern statistics and machine learning. In many applications, data naturally concentrate around lower-dimensional manifolds, such as curves and surfaces, driving the need for methods of non-linear dimension reduction [33, 35] and manifold fitting [14, 36, 37]. Early approaches to modeling these structures, such as principal curves and surfaces [19], provided frameworks by seeking manifolds that pass through the middle of a data cloud to minimize projection error. While subsequent advancements theoretically guaranteed the existence of these curves by imposing global regularization constraints, such as bounding the total length [20] or curvature [34], these formulations remain tied to global optimization heuristics. Because they rely on minimizing projection errors subject to rigid structural penalties, they often struggle to flexibly adapt to variations in local data density. As a compelling alternative, the concept of the density ridge emerged to define the hidden skeleton of a dataset directly through the local differential geometry of the underlying probability density 1

function. By tracking the regions of highest density constrained by local subspaces, density ridges form low-dimensional structures representing the backbone of the data. Rather than imposing arbitrary global length or curvature constraints, they capture the local geometry determined by the data’s distribution. In practice, the Subspace Constrained Mean Shift (SCMS) algorithm, developed by [26], was designed to extract these ridges from data. Prominent applications include the identification of cosmic web filaments and large-scale structures in astrophysics [8], the extraction of seismic fault lines from earthquake data [23], the tracing of anatomical structures such as blood vessels in medical imaging [11], and the optimization of geospatial police patrol routes by mapping crime hot-spot filaments [25]. In this paper, we study the theoretical properties of the SCMS algorithm by clarifying the definition of the ridges it targets, and we quantify its performance in terms of both statistical consistency and computational complexity.

1.1

The Static Ridge and the SCMS Algorithm

We first introduce the concept of density ridges used in the literature. Consider a probability density function f ∈ C 2 (Rd , R) for d ≥ 2. Let k ∈ {1, . . . , d−1} be a fixed integer representing the dimension of the target manifold. For any point x ∈ Rd , let λ1 (x) ≥ · · · ≥ λd (x) be the sorted eigenvalues of the Hessian ∇2 f (x). Let E∥ (x) denote the k-dimensional subspace spanned by the eigenvectors corresponding to the top k eigenvalues of ∇2 f (x), and let E⊥ (x) denote the (d − k)-dimensional subspace spanned by the trailing d − k eigenvectors. We refer to E∥ (x) as the leading eigenspace and E⊥ (x) as the trailing eigenspace of the density Hessian. Let V (x) be the d × (d − k) matrix whose columns form an orthonormal basis for the trailing eigenspace E⊥ (x). We define the orthogonal projection matrix onto this trailing eigenspace as Π(x) = V (x)V (x)⊤ , and the projected density gradient onto this trailing eigenspace is denoted ξ(x) = Π(x)∇f (x). The ridge concept uses the gradient and Hessian (and its eigenvalues and eigenvectors) of the density, and traces its origins to the image processing and computer vision literature [11]. We will refer to such a ridge concept as the static ridge, defined as  Rstatic (f ) = x ∈ Rd : ξ(x) = 0 and λk+1 (x) < 0 . (1) While this is often referred to in the literature as k-dimensional ridges (or k-ridges), we treat k as a fixed parameter throughout this paper and omit it from our notation for brevity. A point belongs to the static ridge if the density gradient is orthogonal to the trailing eigenspace, where the density is a local maximum along any direction in the trailing eigenspace. This geometric interpretation can be seen from the fact that the eigenvalues represent the second directional derivatives of the density along fixed eigenvector directions. Hence the static ridge extends the concept of a local maximum to a set of local maxima constrained in particular subspaces, which forms a k-dimensional structure. In practice, the population density f is unknown and should be estimated from a finite sample X1 , . . . , Xn ∼ f . This is typically achieved using a Kernel Density Estimator (KDE) with a smoothing kernel K and bandwidth h > 0. In the statistical literature, there are a series of papers on the statistical inference for the static ridge using the plug-in estimator based on KDE [6, 16, 28–30]. To extract the static ridge from the estimated density, the SCMS algorithm developed in [26] also utilizes the geometric intuition above for the static ridge. The algorithm builds upon the seminal Mean Shift algorithm, which is a nonparametric mode-seeking technique widely used in clustering and computer vision [9, 10, 15]. The Mean Shift operator, denoted MSn,h : Rd → Rd , at each point x computes a vector pointing toward the locally weighted mass center using neighboring data points of x. The statistical properties of Mean Shift algorithm and relevant hill-climbing algorithms for clustering have been studied in [1–5, 7, 32].

2

In particular, [26] utilize a Gaussian kernel, under which the Mean Shift vector becomes proportional to the gradient of the estimated log-density, that is, MSn,h (x) ∝ ∇ log fb(x). The SCMS algorithm iteratively updates a test point xm by projecting this Mean Shift vector onto an estimated trailing eigenspace. In their original formulation, this trailing eigenspace is determined by the eigenvectors of the log-density Hessian, ∇2 log fb(x). Denoting the corresponding projection matrix as b log (x), the update rule for the SCMS algorithm in [26] is Π b log (xm )MSn,h (xm ). xm+1 = xm + Π

(2)

See Figure 1 for an illustration of the SCMS algorithm’s performance in ridge estimation.

Figure 1: Visualization of the static ridge and SCMS. Left panel: A 3D perspective plot of the population density f , modeled as a 3-component Gaussian mixture, where the solid red line represents the static ridge. Right panel: SCMS trajectories. The grey dots are data samples drawn from the population model, and the blue contours represent the resulting KDE. Green dots indicate the starting points for the SCMS algorithm. The orange points illustrate the path of the algorithm. The trajectories are seen to converge onto the estimated ridge of the log-KDE (in red). Because both the Mean Shift operator and the projection matrix are based on log fb(x), the SCMS algorithm aims for the ridge based on the log-density rather than the raw density as the target. Despite this distinction, it became a widely held belief in the literature [16, 17, 38, 39] that the trajectories of SCMS converged to the static ridge Rstatic (f ) (or its log-density equivalent) in the asymptotic setting when n → ∞ and m → ∞. However, a recent discovery by [31] demonstrated that the relationship between SCMS and the static ridge is more complex than previously assumed. As illustrated by their counterexample (also presented in Appendix A.1.1), the SCMS algorithm does not necessarily converge to the static ridge. The intuition behind this divergence lies in the nature of the SCMS vector field: as a point flows along the trajectory of ξ, the trailing eigenspace Π itself dynamically rotates. Because the static ridge is defined by pointwise directional derivatives, it fails to account for the continuous rotation of the projection matrix along a flow line. To resolve this, [31] developed new algorithms that have a statistical guarantee to converge to the static ridge Rstatic (f ). This still leaves a critical open problem: if the SCMS algorithm does not estimate the static ridge, what structure is it actually finding?

3

1.2

The Stable Ridge: From Geometry to Dynamics

In this paper, we propose a paradigm shift by introducing the concept of the stable ridge, and proving that it is this structure that the SCMS algorithm actually targets and extracts. Unlike the static ridge, which is defined by the Hessian ∇2 f (x), the stable ridge is defined by the topology of the dynamical system governed by ξ. To ensure the stability of this system is welldefined, we assume the density f ∈ C 3 (Rd , R). This higher-order smoothness is necessary because the Jacobian Jξ := ∇ξ incorporates the derivatives of the projection matrix Π. For any x ∈ Rd , let µ1 (x), . . . , µd (x) denote the eigenvalues of this (asymmetric) Jacobian sorted by their real parts in descending order, that is, Re(µ1 (x)) ≥ · · · ≥ Re(µd (x)), where for a complex number z, Re(z) denotes its real part. We define the stable ridge as  Rstable (f ) = x ∈ Rd : ξ(x) = 0 and Re(µk+1 (x)) < 0 . (3) Compared to the static ridge definition in (1), the stable ridge differs only by replacing the condition λk+1 (x) < 0 with Re(µk+1 (x)) < 0. This subtle spectral shift, however, reflects a fundamental conceptual transition from a static geometric perspective to a dynamical systems approach. In classical optimization, local attractors are defined by the zeros of the gradient vector field (i.e., ∇f (x) = 0), with stability governed by the negative eigenvalues of its Jacobian, which is the Hessian ∇2 f (x). Our dynamic framework mirrors this structure for the vector field defined by the projected gradient ξ: the stable ridge is defined by ξ(x) = 0, with stability governed by the negative real parts of the eigenvalues of its Jacobian, Jξ (x). This extension can be clearly illustrated by examining the extreme case where the dimension is k = 0. Under the static definition (1), setting k = 0 requires all d eigenvalues of the Hessian to be negative (λ1 (x) < 0), which recovers the standard second-order condition for a local maximum (or mode). Under our stable definition (3), k = 0 implies the trailing eigenspace spans the entirety of Rd , such that the projection matrix Π(x) becomes the identity matrix Id , the projected gradient ξ(x) reduces to the raw gradient ∇f (x), and Jξ (x) is simply the Hessian ∇2 f (x). Thus, the stable ridge also collapses to the set of local maxima. However, this comparison highlights the fundamental shift in perspective: while the static definition views the mode simply as a peak, the stable definition characterizes it as an equilibrium point that acts as an attractor from all d directions in the dynamical system. However, when the dimension k ≥ 1, these two concepts of static and stable ridges are generally not equivalent. Because the SCMS algorithm can be viewed as a discretization of a continuous dynamical system, the stable ridge—rather than the static ridge—serves as the correct theoretical target for its convergence, the rigorous proof of which is the overarching goal of this paper. A detailed mathematical derivation quantifying the divergence and coincidence between the static and stable ridge definitions (specified when d = 2 and k = 1) with illustrative examples, is provided in Appendix A. Paper Organization. The remainder of this paper is organized as follows. Section 2 introduces the ridge-regular function class as part of our core assumptions for the model. Section 3 characterizes the geometric properties of the stable ridge, studies the continuous flow driven by the vector field ξ, and shows its capability of surjective recovery of the full stable ridge. Section 4 transitions to the discrete algorithmic setting, where we develop the generalized SCMS algorithm (still based on the true density) using a constant step size, and establish its uniform R-linear convergence and topological surjectivity. Section 5 introduces the sample version of the generalized SCMS based on KDE and derives the convergence rate of the Hausdorff error in estimating the true stable ridge. Section 6 bridges our generalized framework back to the original SCMS algorithm operating on the log-density. There, we expose the polynomial-time computational complexity of the original SCMS, which is caused by implicitly coupling the step size to the smoothing bandwidth via the Mean Shift 4

operator, and demonstrate how our generalized framework with a constant step size provides a more efficient solution. Section 7 provides numerical results validating both the statistical consistency and the computational complexities derived in the text. Finally, Appendix A provides examples demonstrating the relationship between the static and stable ridge definitions, and Appendix B contains the proofs for all the theoretical results in this paper.

2

Assumptions and the Ridge-Regular Class

Let D ⊂ Rd be a compact and connected domain of interest for ridge extraction, so that we focus our analysis on the restricted stable ridge R(f ) defined as R(f ) = Rstable (f ) ∩ D = {x ∈ D : ξ(x) = 0 and Re(µk+1 (x)) < 0} .

(4)

We introduce some regularity assumptions used in this paper. Assumption (A1). The density function satisfies f ∈ C 4 (Rd , R) with its partial derivatives up to the fourth order uniformly bounded over Rd . Furthermore, there exists a constant c0 > 0 such that f (x) ≥ c0 for all x ∈ D. While C 3 smoothness is sufficient to define the stable ridge, we require C 4 smoothness to study the consistency and convergence of the SCMS algorithms. Note that while smoothness is formulated on Rd to accommodate infinite-support kernels such as the Gaussian kernel, if we choose a kernel with compact support, then this assumption only needs to hold within a neighborhood of D. To define the geometric conditions required for stable ridge extraction, we collect the requirements into a function class. g For any function g ∈ C 3 (Rd , R), let the quantities, sets and operators λgj (x), E⊥ (x), Πg (x), ξg (x), g and µj (x) be defined analogously to those of f . For any ϵ ≥ 0, define Rϵ (g) = {x ∈ D : ∥ξg (x)∥ ≤ ϵ and Re(µgk+1 (x)) < 0}, such that R0 (f ) = R(f ). Definition 1. The function g is said to belong to the ridge-regular class F (η, γ, κ, D) for some constants η, γ, κ > 0, if it satisfies the following three properties for all x ∈ D. 1. λgk (x) − λgk+1 (x) ≥ η. 2. The set Rκ (g) is contained within the interior of D denoted by int(D), and g sup{u⊤ Jξg (x)u : u ∈ E⊥ (x), ∥u∥ = 1} ≤ −γ.

 3. inf ∥ξg (x)∥ : x ∈ D, Re(µgk+1 (x)) ≥ − γ2

≥ κ.

We assume the population density satisfies these geometric constraints. Assumption (A2). f ∈ F(η, γ, κ, D) for some positive constants η, γ, and κ. Remark 1. The properties required for f ∈ F (η, γ, κ, D) guarantee the statistical and computational stability of the SCMS algorithms. Property 1 provides the necessary spectral gap to ensure the projection matrix Π is uniquely defined and smooth on D. Property 2 forces the symmetric part of the Jacobian to be negative definite when restricted to the trailing eigenspace; this guarantees the monotonic decay of ∥ξ∥ along the SCMS trajectories. Finally, Property 3 regularizes the vector field ξ far away from the ridge to prevent the emergence of spurious features under small perturbations. 5

Notation and Global Constants. To streamline the theoretical analysis in the subsequent sections, we define the following global constants. L = sup ∥Jξ (x)∥

and

M=

x∈D

1 sup ∥∇2 ξ(x)∥, 2 x∈D

where for a matrix A, ∥A∥ denotes its induced spectral norm; for a vector field v : Rd → Rd , its Hessian ∇2 v(x) is a bilinear map, and we define its multilinear operator norm as ∥∇2 v(x)∥ = sup∥u1 ∥=1,∥u2 ∥=1 ∥∇2 v(x)[u1 , u2 ]∥. Based on Assumptions (A1) and (A2), it is straightforward to see L < ∞ and M < ∞. Furthermore, we define d∂ = inf {∥x − y∥ : x ∈ Rκ (f ), y ∈ ∂D} . It can be shown that Rκ (f ) is a compact set (see the proof of Lemma 1 below). Hence d∂ is positive based on Property 2 of F (η, γ, κ, D). We collect some notation used throughout the paper. Let A and B be two non-empty subsets of Rd . The Hausdorff distance between A and B is dH (A, B) = max {d(A|B), d(B|A)} ,

(5)

where d(A|B) = supx∈A inf y∈B ∥x − y∥ and d(B|A) = supy∈B inf x∈A ∥x − y∥. For any x ∈ Rd and r > 0, let Br (x) denote the closed Euclidean ball of radius r centered at x, that is, Br (x) = {y ∈ Rd : ∥y − x∥ ≤ r}. Finally, for positive sequences an and bn , we write an = O(bn ) if there exists a constant C > 0 such that an ≤ Cbn for all n large enough, and we write an ≍ bn if both an = O(bn ) and bn = O(an ) hold.

3

Geometry of the Stable Ridge and Continuous SCMS

In this section, we describe the geometry of the true stable ridge R(f ) and analyze the behavior of the continuous flow determined by the vector field ξ. This continuous flow serves as the infinitesimal limit of the SCMS algorithm at the population level. By understanding the properties of this dynamical system, we lay the groundwork for analyzing the discretization of the flow and the algorithm’s convergence and statistical consistency in subsequent sections. Before tracking the flow, we first confirm that the stable ridge as our target structure is a geometrically well-behaved object under the ridge-regular assumptions. Lemma 1. Under Assumptions (A1) and (A2), the stable ridge R(f ) is a compact, C 2 -smooth, k-dimensional submanifold of Rd without boundary. With the regularity of the stable ridge established, we turn our attention to the trajectories of the continuous population version of SCMS. Let φt (x) denote the flow driven by the vector field ξ, with the initial value x, that is, φ̇t (x) :=

d φt (x) = ξ(φt (x)), dt

We define its total arc length as

with φ0 (x) = x.

(6)

Z ∞ L(x) =

∥ξ(φt (x))∥ dt.

(7)

0

As a practical ridge estimation technique, it is crucial for the paths of continuous SCMS initialized near the target ridge remain nearby and converge to the ridge within a finite distance. The following lemma demonstrates that the geometric conditions of the ridge-regular class guarantee an exponential decay in the magnitude of the vector field along the paths, ensuring this convergent behavior. 6

Lemma 2. Suppose Assumptions (A1) and (A2) hold. For ϵ ∈ (0, κ] and x ∈ Rϵ (f ), the flow φt satisfies ∥ξ(φt (x))∥ ≤ e−γt ∥ξ(x)∥ ≤ ∥ξ(x)∥ ≤ ϵ for all t ≥ 0. (8) Consequently, φt (x) ∈ Rϵ (f ) for all t ≥ 0 and L(x) ≤

1 ∥ξ(x)∥. γ

(9)

Furthermore, φt (x) converges to a limit point in R(f ) as t → ∞, which we denote Φ(x), i.e., Φ(x) = lim φt (x), for x ∈ Rϵ (f ). t→∞

(10)

Because any point x ∈ Rϵ (f ) flows onto the ridge R(f ), it is informative to quantify how tightly the neighborhood Rϵ (f ) bounds the true ridge. The following corollary establishes that their Hausdorff distance shrinks linearly with ϵ. The corollary will be useful in quantifying the perturbation of the stable ridge, as a slightly perturbed ridge remains in Rϵ (f ) for some small ϵ > 0.
γ Corollary 1. Suppose Assumptions (A1) and (A2) hold. Let ϵ satisfy 0 ≤ ϵ ≤ min κ, Lγ M , Lκ . Then we have 1 1 ϵ ≤ dH (Rϵ (f ), R(f )) ≤ ϵ. (11) 2L γ While Lemma 2 ensures that every point in the neighborhood converges to the ridge, it does not guarantee that the entire ridge is covered by these limit points. To ensure our extraction process discovers the complete manifold, we need to prove the surjectivity of limit map Φ defined by the flow φt . For ϵ ∈ (0, κ), we define the boundary of the neighborhood as ∂Rϵ (f ) = {x ∈ D : ∥ξ(x)∥ = ϵ and Re(µk+1 (x)) < 0}.

(12)

The following theorem establishes that initializing the flow φt from ∂Rϵ (f ) is sufficient to recover every point on the target ridge R(f ). Theorem 1. Suppose Assumptions (A1) and (A2) hold. For ϵ ∈ (0, κ), the limit map Φ defined in (10) is surjective from ∂Rϵ (f ) onto R(f ). That is, for every point x∗ ∈ R(f ), there exists at least one starting point x0 ∈ ∂Rϵ (f ) such that Φ(x0 ) = x∗ . Finally, for the limit map Φ(x) to be useful in bounding the errors in subsequent discrete and finite-sample settings, it must be stable under small perturbations. We conclude this section by showing the the continuous differentiability of Φ. Lemma 3. Suppose Assumptions (A1) and (A2) hold. The limit map Φ defined in (10) is continuously differentiable on Rκ (f ). Furthermore, there exists a uniform constant LΦ > 0 such that for all x ∈ Rκ (f ), we have ∥JΦ (x)∥ ≤ LΦ , (13) where JΦ (x) is the Jacobian of Φ.

7

4

Stability and Convergence of the Generalized SCMS Algorithm

To bridge the gap between the continuous population version of SCMS studied in the previous section and the practical implementation of the algorithm, we analyze the discretized approximation of the flow φt . To distinguish our mathematical framework from the original formulation by [26], which couples the step size to the smoothing bandwidth in the algorithm and is implicitly based on the log-density, we define the generalized SCMS algorithm. This generalized version operates directly on the vector field ξ(x) derived from the density using a constant step size α > 0 specified by the user. We define the one-step operator Gα : D → Rd for the discrete algorithm as Gα (x) = x + αξ(x).

(14)

The generalized SCMS sequence initiated at x0 is generated by the m-fold composition of this operator, denoted xm = Gm α (x0 ), for m = 1, 2, . . . . By imposing appropriate step-size bounds derived from our global constants that are defined in Section 2, we establish that the discrete map Gα prevents the sequence from escaping the neighborhood of the ridge. Furthermore, we prove that the magnitude of projected gradient ξ decays at a geometric (linear) rate along the path in the algorithm, which recovers the entire stable ridge, mirroring what is achieved in the continuous setting. The following lemma provides the required step-size upper bounds. We define the total arc length of the discrete path generated by the sequence {xm }∞ m=0 as Lα (x0 ) =

∞ X

∥xm+1 − xm ∥.

m=0

Lemma 4. Suppose Assumptions (A1) and (A2) hold. Let ϵ ∈ (0, κ]. If the sequence {xm }∞ m=0 is initialized at x0 ∈ Rϵ (f ) and the step size satisfies   γ d∂ γ , , , (15) α < min L2 4M ϵ ϵ then for all m ≥ 0, we have (a) xm ∈ Rϵ (f ). (b) ∥ξ(xm+1 )∥ ≤ ρ∥ξ(xm )∥, where ρ = 1 − αγ 4 < 1. (c) Lα (x0 ) ≤ γ4 ∥ξ(x0 )∥. With the length of the discrete trajectory bounded and the projected gradient decaying at a geometric rate, we can now show that the generalized SCMS algorithm converges to the stable ridge at a uniform rate. Theorem 2. Suppose Assumptions (A1) and (A2) hold, and the step-size condition (15) of Lemma 4 is satisfied. For any ϵ ∈ (0, κ] and any initial point x0 ∈ Rϵ (f ), the sequence {xm }∞ m=0 converges to a limit point x∞ ∈ R(f ). Furthermore, the convergence is R-linear uniformly for all starting points x0 ∈ Rϵ (f ), that is, for all m ≥ 0, ∥xm − x∞ ∥ ≤ Cρm , αγ where C = 4ϵ γ and ρ = 1 − 4 .

8

(16)

Theorem 2 guarantees the existence of a well-defined map, which we denote Φα : Rϵ (f ) → R(f ), given by the limit of the iterative sequence: Φα (x) = lim Gm α (x), for x ∈ Rκ (f ). m→∞

(17)

To effectively utilize this map in the subsequent finite-sample analysis, we first show that it is continuous. Proposition 1. Suppose Assumptions (A1) and (A2) hold, and the step-size condition (15) of Lemma 4 is satisfied. For any ϵ ∈ (0, κ], the limit map Φα (x) is continuous on Rϵ (f ). Furthermore, it is important to quantify how closely this map Φα based on discrete trajectories approximates Φ based on continuous flows and defined in (10). The next proposition demonstrates that their difference is linearly bounded by the step size α. Proposition 2. Suppose Assumptions (A1) and (A2) hold. Let ϵ ∈ (0, κ) and x0 ∈ Rϵ (f ). Assume the step size α satisfies the condition (15) in Lemma 4 and additionally satisfies α<

κ−ϵ . Lϵ

Then we have

 ∥Φα (x0 ) − Φ(x0 )∥ ≤

 2LΦ L ∥ξ(x0 )∥ α. γ

(18)

(19)

Finally, mirroring the surjectivity for the continuous flow established in Theorem 1, we can prove that the generalized SCMS algorithm is also capable of recovering the entire stable ridge, when initialized from the outer boundary of the neighborhood. Theorem 3. Suppose Assumptions (A1) and (A2) hold. For ϵ ∈ (0, κ], let the step size α > 0 satisfy   γ γ 1 d∂ α < min , , , . (20) L2 4M ϵ L ϵ Then the limit map Φα is surjective from ∂Rϵ (f ) onto R(f ). That is, for every point x∗ ∈ R(f ), there exists a starting point x0 ∈ ∂Rϵ (f ) such that Φα (x0 ) = x∗ .

5

Statistical Consistency of the Empirical Algorithm

Having characterized the geometric and topological stability of the discrete algorithm based on the population density, we now turn to the practical setting where the true density f is unknown. In practice, given a finite set of independent and identically distributed sample points X1 , . . . , Xn ∼ f , we can estimate f by using the Kernel Density Estimator (KDE) over Rd :   n 1 X x − Xi b f (x) = K , nhd i=1 h where K : Rd → R is the smoothing kernel and h ≡ hn > 0 is the bandwidth. Using this estimator, we can estimate the gradient ∇f (x) and Hessian ∇2 f (x) by ∇fb(x) and 2b b1 (x) ≥ · · · ≥ λ bd (x) denote the sorted eigenvalues of ∇2 fb(x). Let E b⊥ (x) ∇ f (x), respectively. Let λ denote the eigenspace spanned by the eigenvectors corresponding to the smallest d − k eigenvalues 9

b of ∇2 fb(x), and let Π(x) denote the orthogonal projection matrix onto this subspace. Then ξ(x) can be estimated by the plug-in estimator b = Π(x)∇ b ξ(x) fb(x). b Furthermore, we define Jξb(x) = ∇ξ(x) as the Jacobian of ξb at x, with its eigenvalues denoted by µ b1 (x), . . . , µ bd (x), sorted by their real parts in a descending order. The plug-in stable ridge estimator is n o b = 0 and Re(b R(fb) = x ∈ D : ξ(x) µk+1 (x)) < 0 . For a given threshold ϵ > 0, corresponding to Rϵ (f ), we define n o b Rϵ (fb) = x ∈ D : ∥ξ(x)∥ ≤ ϵ and Re(b µk+1 (x)) < 0 . b that is, Similarly, let φ bt (x) denote the trajectory driven by the vector field ξ, d bφ φ bt (x) = ξ( bt (x)), dt

with φ b0 (x) = x.

The empirical one-step SCMS operator, parameterized by a constant step size α > 0, is defined as b b α (x) = x + αξ(x). G The sequence for the sample version of the generalized SCMS algorithm is generated by the m-fold bm x0 ). composition x bm = G α (b In this section, we derive the statistical consistency and convergence of this empirical algorithm, where the rates of convergence in the error bounds are quantified in terms of the sample size n, bandwidth h, and step size α.

5.1

Assumptions for Statistical Estimation

To guarantee the uniform consistency of fb and its higher-order derivatives, we require some regularity conditions on the kernel and bandwidth. For a multi-index a = (a1 , . . . , ad ) of non-negative integers Pd |a| f with |a| = i=1 ai = m, let Da f (x) = ∂xa∂1 ...∂x ad (x) denote the corresponding partial derivative of 1 d order m. We require the following standard conditions for kernel density estimation and its derivatives. Assumption (B1). The kernel K ∈ C 4 (Rd , R) is a symmetric probability density function with R bounded partial derivatives up to the fourth order. It satisfies uu⊤ K(u)du = κ2 I, where 0 < κ2 < ∞. Furthermore,     x−y a d : h > 0, y ∈ R , |a| = m for m = 0, . . . , 4, Fm = x 7→ D K h form Vapnik-Chervonenkis (VC) subgraph classes. Assumption (B2). As the sample size n → ∞, the bandwidth sequence hn satisfies hn → 0 and nhd+8 n → ∞. log n 10

Remark 2. Assumptions (B1) and (B2) are mild for nonparametric estimation of densities and their higher-order derivatives. In (B1), the C 4 smoothness guarantees that the vector field ξb is sufficiently smooth on D. The requirement that the families Fm form VC subgraph classes is standard for the L∞ stochastic error bounds in empirical process theory [12, 18]. Assumption (B2) places a lower bound on the bandwidth decay rate specifically to ensure that the uniform stochastic error of the fourth-order derivative estimator vanishes uniformly over the domain D. The following lemma provides the uniform convergence rates for the KDE and its derivatives up to the fourth order, which serve as the baseline for all subsequent analyses for the plug-in estimators based on KDE and its derivatives. We restrict our statistical analysis to the compact domain D given in Assumption (A1). Lemma 5. Under Assumptions (A1), (B1), and (B2), for any constant c > 0, there exist positive constants C0 , C1 , C2 , C3 , C4 and a deterministic function ρ4 (h) = o(1) such that for all n sufficiently large, the following bounds hold simultaneously with probability at least 1 − n−c . ! r log n ab a 2 =: C0 r0 (n), (21) max sup |D f (x) − D f (x)| ≤ C0 h + nhd |a|=0 x∈D ! r log n ab a 2 max sup |D f (x) − D f (x)| ≤ C1 h + =: C1 r1 (n), (22) nhd+2 |a|=1 x∈D ! r log n ab a 2 =: C2 r2 (n), (23) max sup |D f (x) − D f (x)| ≤ C2 h + nhd+4 |a|=2 x∈D ! r log n max sup |Da fb(x) − Da f (x)| ≤ C3 h + =: C3 r3 (n), (24) nhd+6 |a|=3 x∈D ! r log n ab a max sup |D f (x) − D f (x)| ≤ C4 ρ4 (h) + =: C4 r4 (n). (25) nhd+8 |a|=4 x∈D

5.2

Ridge Regularity for the Estimated Density

Recall that the SCMS algorithm operates dynamically based on the density gradient and the eigenspaces of the density Hessian. We need to translate the perturbation bounds for the density derivatives in Lemma 5 into perturbation bounds for the geometric quantities defining SCMS. To ensure the sample version of the generalized SCMS algorithm behaves the same as its population version, for any constant c > 0, we define En (c) as the event where the uniform convergence bounds in Lemma 5 hold simultaneously. By Lemma 5, it is guaranteed that P(En (c)) ≥ 1 − n−c for all sufficiently large n. Lemma 6. Under Assumptions (A1)–(A2) and (B1)–(B2), for any c > 0, there exists an integer N0 such that for all n ≥ N0 , conditionally on the event En (c), the following bounds hold. bj (x) − λj (x)| ≤ dC2 r2 (n). (a) supx∈D maxj |λ b (b) supx∈D ∥Π(x) − Π(x)∥ ≤ 2dη −1 C2 r2 (n). (c) There exists a constant Cξ > 0 such that b − ξ(x)∥ ≤ ϵξ (n) := Cξ r2 (n). sup ∥ξ(x) x∈D

11

(d) There exists a constant CJ > 0 such that sup ∥Jξb(x) − Jξ (x)∥ ≤ ϵJ (n) := CJ r3 (n).

x∈D

Consequently, for a sufficiently large sample size, the geometric quantities defining SCMS based on fb are stochastically close to those based on the true density f . The following corollary formalizes that fb satisfies the ridge-regular condition for f as defined in Definition 1, with parameters scaled by a factor of 1/2 to absorb the uniform stochastic errors. Corollary 2. Assume Assumptions (A1)–(A2) and (B1)–(B2) hold with f ∈ F(η, γ, κ, D). For any c > 0, there exists an integer N1 ≥ N0 such that for all n ≥ N1 , conditionally on En (c), fb satisfies Assumption (A1) and η γ κ  , , ,D . fb ∈ F 2 2 2 Specifically, this guarantees the following. (a) fb ∈ C 4 (Rd , R), and there exist constants Bℓ > 0 such that max|a|=ℓ supx∈D |Da fb(x)| ≤ Bℓ , for ℓ = 0, . . . , 4. bk (x) − λ bk+1 (x) ≥ η > 0 for all x ∈ D. (b) λ 2 (c) For all x ∈ Rκ/2 (fb), b⊥ (x), ∥u∥ = 1} ≤ − sup{u⊤ Jξb(x)u : u ∈ E  b (d) inf ∥ξ(x)∥ : x ∈ D, Re(b µk+1 (x)) ≥ − γ4

5.3

γ < 0. 2

≥ κ2 > 0.

Consistency of the Plug-in Ridge Estimator

b α , which are expected to converge to the plug-in Before analyzing discrete trajectories induced by G b ridge estimator R(f ), we first confirm that the target R(fb) is close to the true ridge R(f ). The following theorem demonstrates that the Hausdorff distance between the estimated ridge and the true ridge shares the same rate of convergence as the estimation of the Hessian. Theorem 4. Assume Assumptions (A1)–(A2) and (B1)–(B2) hold. For any c > 0, there exists an integer N2 ≥ N1 such that for all n ≥ N2 , conditionally on the event En (c), ! r 2ϵ (n) log n ξ 2 dH (R(fb), R(f )) ≤ =O h + . (26) γ nhd+4 Remark 3. Minimizing the Hausdorff error rate in Theorem 4 requires  balancing the squared bias p rate O(h2 ) with the uniform stochastic error rate O log n/(nhd+4 ) . This produces the balancing 1

d+8

bandwidth h ≍ (log n/n) d+8 such that nh log n ≍ 1, which however narrowly fails the requirement in Assumption (B2) for the uniform convergence of the fourth derivatives. To address this, we apply an arbitrarily small logarithmic correction. For any constant δ > 0, choosing the bandwidth  h≍

(log n)1+δ n 12

1  d+8

d+8

δ ensures that nh log n ≍ (log n) → ∞, which also satisfies Assumption (B2). Under this choice, the estimator is slightly oversmoothed, and the Hausdorff error is dominated by the bias, which yields   2 ! (log n)1+δ d+8 b . dH (R(f ), R(f )) = O n

5.4

Convergence of the Generalized SCMS Algorithm: Sample Version

We now focus on the extraction of the estimated ridge R(fb) using the empirical SCMS operator b b α (x) = x + αξ(x). G For any ϵ ∈ (0, κ/2], we define the boundary of Rϵ (fb) as  b ∂Rϵ (fb) = x ∈ D : ∥ξ(x)∥ = ϵ and Re(b µk+1 (x)) < 0 . (27) b bm b Let Xbm = {G α (x) : x ∈ ∂Rϵ (f )} denote the image of this boundary after m iterations of Gα . The following theorem mirrors the population result in Theorem 2, proving that the sample version of the generalized SCMS algorithm uniformly converges to the estimated ridge at an Rlinear rate. Theorem 5. Let ϵ ∈ (0, κ/2]. Under Assumptions (A1)–(A2) and (B1)–(B2), for any c > 0, there exists an integer N2 ≥ N1 such that for all n ≥ N2 , conditionally on the event En (c), the following holds. For any point x b0 ∈ ∂Rϵ (fb), if the step size α > 0 satisfies   γ γ d∂ α < min , , , (28) 8L2 16M ϵ ϵ bm x0 )}∞ b∞ ∈ R(fb), and then the sequence {b xm = G α (b m=0 converges to a limit point x 8ϵ  αγ m ∥b xm − x b∞ ∥ ≤ 1− , for all m ≥ 0. γ 8

(29)

Just as in the population setting, it is insufficient to characterize the behavior of the sample version of the generalized SCMS algorithm by pointwise convergence alone, as it is possible that the algorithm only recovers an arbitrary fraction of the stable ridge. We need to demonstrate topological surjectivity onto the estimated stable ridge to ensure a full recovery of the target. Theorem 6. Let ϵ ∈ (0, κ/2]. Under Assumptions (A1)–(A2) and (B1)–(B2), for any c > 0, there exists an integer N2 ≥ N1 such that for all n ≥ N2 , conditionally on the event En (c), the following holds. If the step size α > 0 satisfies   γ γ 1 d∂ , α < min , , , (30) 8L2 16M ϵ 2L ϵ b b α (x) = limm→∞ G bm then the limit map Φ α (x) is surjective from the boundary set ∂Rϵ (f ) onto the estimated ridge R(fb). Corollary 3. Assume the conditions of Theorem 6 hold. For any c > 0, there exists an integer N2 such that for all n ≥ N2 , with probability at least 1 − n−c , we have 8ϵ  αγ m dH (Xbm , R(fb)) ≤ 1− . (31) γ 8 13

We are now ready to state the capstone result of this section. By combining Theorem 4 and Corollary 3 and using a triangle inequality for the Hausdorff distance, we provide the total error bound for ridge estimation using the sample version of the generalized SCMS algorithm. Theorem 7. Assume the conditions of Theorem 6 hold. For any c > 0, there exists an integer N2 such that for all n ≥ N2 , with probability at least 1 − n−c , we have αγ m 2ϵξ (n) 8ϵ  1− dH (Xbm , R(f )) ≤ + . γ 8 γ

(32)

Remark 4. In Theorem 7, the total error in estimating the population stable ridge is decomposed into the computational error and the statistical error, corresponding to the first and second terms on the right side of (32), respectively. The computational error vanishes exponentially fast with the number of iterations m, while the statistical error depends on the sample size n and bandwidth h. To achieve the convergence rate for ridge estimation carefully balanced in Remark 3, we can b α until the computational error is dominated by the statistical error. Using execute the operator G  1 2  (log n)1+δ
d+8 (log n)1+δ
d+8 , the statistical error is bounded by ϵ (n) = O , the bandwidth h ≍ ξ n n following from Lemma 6. m By setting the computational error 8ϵ ≍ ϵξ (n), we find the required number of γ (1 − αγ/8) iterations before stopping the algorithm is      log(8ϵ/γ) − log(ϵξ (n)) 1 m∗ = = O log = O(log n). (33) − log(1 − αγ/8) ϵξ (n) Evaluating at this stopping time, we obtain the rate of convergence for the total error using the sample version of the generalized SCMS algorithm as  2 !  (log n)1+δ d+8 b dH (Xm∗ , R(f )) = O . n

6

Convergence and Computational Complexity of the Original SCMS Algorithm

The generalized SCMS framework developed in Sections 4 and 5 operates on the density f and its estimator fb using a constant step size α. We now extend the theory there to the original SCMS algorithm [26], which is based on an estimator of the log-density and couples its step size to the smoothing bandwidth. Define the log-density as p(x) = log f (x) and its plug-in estimator as pb(x) = log fb(x). While 1 the gradients of the density and log-density have the same direction (∇p(x) = f (x) ∇f (x)), their Hessians have different eigenspaces in general: ∇2 p(x) =

∇2 f (x) ∇f (x)∇f (x)⊤ − . f (x) f (x)2

(34)

This makes the log-density ridge R(p) sometimes distinct from the density ridge R(f ). As discussed in Section 1, the theoretical target of the original SCMS algorithm is the log-density ridge R(p). To establish the convergence of the original SCMS algorithm, we update Assumption (A2) to Assumption (A2′ ). p ∈ F(ηlog , γlog , κlog , D) for some positive constants ηlog , γlog , and κlog . 14

Furthermore, in (B1), we additionally assume that the smoothing kernel K is radially symmetric. Specifically, we assume there exists a differentiable profile function k : [0, ∞) → [0, ∞) such that the kernel can be written as K(x) = ck k(∥x∥2 ), where ck > 0 is a normalization constant such that K integrates to 1 over Rd . We require that k(u) is non-increasing for all u ≥ 0, and strictly decreasing for all u such that k(u) > 0. The updated Assumption (B1) is called Assumption (B1′ ). We retain Assumptions (A1) and (B2). Because the profile function is non-increasing, its negative derivative yields a non-negative weight profile h(u) = −k′ (u) ≥ 0, which is positive on the support of the kernel (i.e., wherever k(u) > 0). We similarly define ch > 0 as the corresponding normalization constant for h. Let Πlog (x) denote the orthogonal projection matrix onto the eigenspace associated with the b log (x) denote its plug-in estimator smallest d−k eigenvalues of the log-density Hessian ∇2 p(x). Let Π 2 2 log log when ∇ p(x) is replaced by ∇ pb(x). Denote ξ (x) = Π (x)∇p(x), and its Jacobian as Jξlog (x) = b log (x)∇b ∇ξ log (x). Similarly, denote ξblog (x) = Π p(x) and let Jξblog (x) be its Jacobian. To formally introduce the original SCMS algorithm, we first define the Mean Shift vector
Pn x−Xi 2 i=1 Xi h h (35) MSn,h (x) = Pn
− x, x−Xi 2 i=1 h h which computes a locally weighted shift pointing toward the center of mass of the neighboring data points around x, using the profile h. The original SCMS algorithm generates a sequence of iterates b log MSn,h , as specified in (2). {xm }∞ m=0 via the update rule determined by Π In what follows, we show that the original SCMS algorithm can be represented in the form of xm+1 = xm + αn (xm )ξblog (xm ),

(36)

where αn is a varying step size depending on the bandwidth h.
x−Xi 2 ch Pn as another kernel density estimator using the profile By defining qb(x) = nh d i=1 h h h, we can write ∇fb(x) MSn,h (x) = CK h2 = αn (x)∇b p(x), (37) qb(x) (x) ch where CK = 2c is a constant only depending on the kernel, and αn (x) = CK h2 fqb(x) . This reveals k that the Mean Shift vector is parallel to the log-density gradient estimator ∇b p(x), scaled by an adaptive step size αn (x). Notably, when utilizing the standard Gaussian kernel, which is the primary focus of the original SCMS [26], the profiles satisfy k(u) = e−u/2 and h(u) = 21 e−u/2 . At the same time, ch = 2ck (yielding CK = 1), and the density estimators are identical (b q (x) ≡ fb(x)). In this special case, the step size αn (x) becomes independent of x, that is, αn (x) ≡ h2 . Plugging (37) into (2), we immediately recover the form of the original SCMS in (36), or equivb log alently, xm+1 = G αn (xm ), where b

blog (x). b log G αn (x) = x + αn (x)ξ

(38)

b α for the sample version of the This one-step operator for the original SCMS is analogous to G generalized SCMS, and analyses are also similar to those in Section 5. First, similar to Lemma 5, under Assumptions (A1), (B1′ ) and (B2), for any constant c > 0, we have P(En′ (c)) ≥ 1 − n−c for n large enough, where En′ (c) is the event such that max sup |Da pb(x) − Da p(x)| ≤ Cℓlog rℓ (n) |a|=ℓ x∈D

15

for ℓ = 0, 1, . . . , 4,

where rℓ (n) are the rates from Lemma 5, and Cℓlog are positive constants. This can be verified by noticing f (x) ≥ c0 > 0 under Assumption (A1), and using a telescoping argument. Consequently, results similar to Lemma 6 and Corollary 2 in the log-density setting hold. In particular, for the perturbation of the vector field, conditionally on En (c) we have sup ∥ξblog (x) − ξ log (x)∥ ≤ ϵlog ξ (n),

x∈D

q   log n 2 . where ϵlog d+4 ξ (n) = O(r2 (n)) = O h + nh b Furthermore, since both f and qb are valid kernel density estimators, they converge uniformly to the population density f , and hence their ratio converges uniformly to 1 across the domain D, where f (x) ≥ c0 > 0. We define the ratio error as fb(x) −1 . b(x) x∈D q

(39)

∆n = sup

Conditionally on En (c), ∆n → 0 as n → ∞, and hence αn (x) ≍ h2 uniformly in D, ensuring the adaptive step size remains stable. For ϵ ∈ (0, κlog /2], we define the boundary set of Rϵ (b p) as  ∂Rϵ (b p) = x ∈ D : ∥ξblog (x)∥ = ϵ and Re(b µlog k+1 (x)) < 0 . The following lemma guarantees that for a sufficiently large sample size n (which also ensures a sufficiently small adaptive step size αn ), the original SCMS algorithm can fully recover the estimated log-density stable ridge R(b p) if starting from the boundary ∂Rϵ (b p) with high probability. Lemma 7. Let ϵ ∈ (0, κlog /2]. Under Assumptions (A1), (A2′ ), (B1′ ), and (B2), for any c > 0, b log there exists an integer N3 such that for all n ≥ N3 , conditionally on the event En (c), the map G αn log log m b b defined in (38) is a local diffeomorphism on Rϵ (b p), and its limit map Φ = limm→∞ (G ) exists and is surjective from the boundary set ∂Rϵ (b p) onto the empirical log-density ridge R(b p).  log b log )m (x) : x ∈ ∂Rϵ (b = (G p) denote the image of the boundary set after running Let Xbm the original SCMS algorithm for m steps. The final theorem below provides the total convergence bounds for the original SCMS algorithm. Theorem 8. Let ϵ ∈ (0, κlog /2]. Under Assumptions (A1), (A2′ ), (B1′ ), and (B2), for any c > 0, there exists an integer N3 such that for all n ≥ N3 , with probability at least 1 − n−c , we have dH (R(b p), R(p)) ≤ and log dH (Xbm , R(b p)) ≤

where ρn = 1 −

2ϵlog ξ (n) γlog

,

8ϵ m ρ , γlog n

(40)

(41)

CK γlog 2 16 h , and consequently, log 8ϵ m 2ϵξ (n) log dH (Xbm , R(p)) ≤ ρn + . γlog γlog

16

(42)

Remark 5. The result in Theorem 8 exposes the severe computational cost of the original SCMS algorithm, for which the step size αn is proportional to h2 , as the sample size n → ∞ and the bandwidth h → 0. Typically, the stopping criterion for the original SCMS algorithm can be set as ∥ξblog (x)∥ < ϵmin log log for a small value of ϵmin > 0. This implies Xbm ⊂ Rϵmin (b p) and hence dH (Xbm , R(b p)) = O(ϵmin ) using an empirical result similar to Corollary 1. By invoking (41) and applying the inequality log(1 − x) ≤ −x for x ∈ (0, 1), we find the required number of iterations is O(h−2 ). To achieve the targeted convergence rate in (40) for nonparametric estimation, we need to iterate ρm is dominated by the statistical the original SCMS algorithm until the computational error γ8ϵ log n q   log n 2 error ϵlog . Setting the two errors equal yields a required number of iterations: ξ (n) = O h + nhd+4 & ∗

m =

log(8ϵ/γlog ) − log(ϵlog ξ (n)) CK h2 γlog /16

' =O

1 log h2

1 ϵlog ξ (n)

!! .

(43)

Because the log-density estimate shares the same uniform convergence rates as the KDE, we substitute 1  d+8  1+δ the bandwidth h ≍ (log n) derived in Remark 3 to satisfy Assumption (B2). Then the required n number of iterations blows up polynomially with the sample size:  2  2(1+δ) m∗ = O n d+8 (log n)1− d+8 . (44) The high computational complexity is mainly due to the Mean Shift operations implicitly using step size αn ≍ h2 . Our generalized SCMS framework developed in Section 4 and Section 5 offers a solution by applying a constant step size α > 0: xm+1 = xm + αξblog (xm ).

(45)

Similar to Remark 4, when the step size α is appropriately calibrated, with the log-density ridge R(p) as the estimation target, the generalized SCMS algorithm in (45) achieves the statistical error rate only requiring O(log n) iterations. This reduced computational complexity makes log-density ridge extraction feasible when the sample size is massive.

7

Numerical Experiments

To validate the theoretical results regarding the computational complexity of SCMS algorithms and their statistical consistency in estimating the log-density stable ridge, we conduct Monte Carlo simulations based on a probability density f , which is constructed by convolving the uniform distribution supported on M and an isotropic Gaussian distribution N (0, σ 2 I2 ) with σ = 0.5, where M is the 1-dimensional circle of radius r = 2 centered at the origin in R2 . The density f is radially symmetric and can be expressed as f (x) = g(∥x∥) for a function g : [0, ∞) → [0, ∞). Such a model is further discussed in Appendix A.1.3. Recall that p denotes the log density log f . For this density, the stable log-density ridge is Rstable (p) = {x ∈ R2 : ∥x∥ = r∗ }, where r∗ ≈ 1.93. For theoretical purposes, we can find a tubular neighborhood D = {x ∈ R2 : ρmin ≤ ∥x∥2 ≤ ρmax } that satisfies Assumption (A2′ ). This requires (log g)′′ < 0 on this region, or roughly 0.42 ⪅ ρmin < r∗ < ρmax < ∞. From this model, we draw i.i.d. random samples Xi = (r cos(Θi ) + Ui , r sin(Θi ) + Vi ) ∈ R2 , for i = 1, . . . , n, where Θi ∼ Uniform(0, 2π), and (Ui , Vi ) ∼ N (0, 0.52 I2 ) with Θi independent of 17

(Ui , Vi ). The sample sizes n range in {103 , 103.5 , 104 , 104.5 , 105 }. To satisfy Assumption (B2) (also see Remark 4), we choose the bandwidth h = C(n−1 log1+δ n)1/(d+8) , with d = 2 and δ = 0.1, where the constant C is calibrated such that h = 0.3 for n = 103 . We use the Gaussian kernel throughout the experiments. Our goal is to estimate R(p) by extracting R(b p). To examine the trajectories of the original and generalized SCMS algorithms, the initial points x0 are selected at (u, 0) ∈ R2 for u ∈ {0.1, 1.5, 3.0}. Specifically, x0 ∈ D when u = 1.5 and u = 3.0, while the initialization at u = 0.1 forces the algorithms to navigate for some distance on the density landscape before entering D. The stopping criterion for all the SCMS algorithms is set to be ∥ξblog (x)∥ < 10−6 . We calculate the required number of iterations for the SCMS algorithms and their statistical errors by averaging over 15 independent replications. As evidenced by the narrow standard error ribbons in Figure 2 and Figure 4 below, this number of replications is sufficient to make stable asymptotic conclusions.

7.1

Computational Complexity and Algorithmic Stability

First, we want to verify the computational complexity of the SCMS algorithms discussed in Remark 5. As illustrated in Figure 2, the iteration count for the original SCMS exhibits an approximate linear trend on the log-log scale. Remark 5 predicts an approximate iteration complexity scaling of O(h−2 ), which is O(n2/(d+8) ) = O(n0.2 ) up to a logarithm factor given d = 2. To only evaluate large sample asymptotic behavior, we fit a linear regression based on the logarithms of the mean number of iterations and the logarithms of the sample sizes n for n ∈ {104 , 104.5 , 105 }. For u ∈ {1.5, 3.0} (corresponding to x0 ∈ D), the fitted log-log slopes are 0.174 and 0.211, respectively, well approximating the theoretical bound. Even for u = 0.1 (corresponding to x0 ∈ / D), the fitted log-log slope is 0.167. In contrast, evaluating the generalized SCMS across a range of step sizes α ∈ {0.1, 0.2, 0.3, 0.4} reveals a computational paradigm shift. For example, focusing on α = 0.2, the fitted log-log slopes for u ∈ {0.1, 1.5, 3.0} are −0.037, −0.032, and −0.038, respectively, which are approximately zero. The slight negative trend even underscores the higher efficiency of the algorithm as the sample size grows. This confirms that decoupling the step size from the bandwidth in the generalized SCMS algorithm bypasses the O(h−2 ) computational bottleneck for the original SCMS algorithm, as discussed in Remark 5. We highlight that similar behaviors are observed between u ∈ {1.5, 3.0} and u = 0.1, just like for the original SCMS algorithm. This indicates that our theoretical analyses requiring assumptions on D can be extended to regions beyond D. See the Discussion section. While the constant step size for the generalized SCMS algorithm has computational benefits compared to the original SCMS algorithm, α needs to be carefully calibrated. As shown in Figure 2, while increasing the step sizes α from 0.1 to 0.3 can speed up the computation for the generalized SCMS, increasing the step size α from 0.3 to 0.4 instead requires more iterations for convergence. An even larger step size could lead to divergence. This is similar to the well-known overshooting problem in gradient ascent algorithms. Therefore, the step size α needs to be carefully tuned, usually by the trial and error method. On the other hand, there is no such a step-size selection problem for the original SCMS algorithm. Figure 3 visualizes the trajectories of the original SCMS algorithm and the generalized SCMS algorithm with α = 0.2 and 0.4, where it can be seen that the original SCMS requires an excessive sequence of small steps before reaching the estimated ridge, while the generalized SCMS with α = 0.4 exhibits a zig-zagging trajectory which overshoots the target during the early stage of the algorithm.

18

Figure 2: Number of iterations to convergence across varying sample sizes n and starting distances u for the original SCMS and Generalized SCMS with step sizes α = 0.1, 0.2, 0.3, 0.4. The plot axes are visualized on a log-log scale.

Figure 3: Trajectories of the original and generalized SCMS algorithms. The black dot denotes the starting point and the solid green line is the estimated stable ridge R(b p). Setting α = 0.4 for the generalized SCMS results in overshooting and underdamped oscillation across the estimated ridge.

19

7.2

Verification of Statistical Consistency

Next, we validate the statistical consistency of the ridge’s plug-in estimator by evaluating the Hausdorff distance between the estimated ridge R(b p) and the population ridge R(p). Here, R(b p) is extracted using the generalized SCMS algorithm with step size α = 0.2. Because we used the stopping criterion ∥ξblog (x)∥ < 10−6 , while the Hausdorff errors are of order between 10−1 and 10−2 (see Figure 4), the numerical errors induced by the SCMS algorithm are dominated by the statistical error and can be safely ignored. As illustrated in Figure 4, the Hausdorff errors exhibit an approximate linear decay on the loglog scale. As presented in Theorem 8 and Remark 5, theoretically statistical error is asymptotically bounded by a rate of O(n−2/(d+8) ) = O(n−0.2 ) up to a logarithmic factor, given d = 2 and our bandwidth selection rule. We fitted a linear regression based on the logarithm of the mean Hausdorff error and the logarithm of the sample sizes across all simulated scenarios, yielding a fitted regression slope of −0.326, which is faster than the theoretical upper bound of −0.2. This occurs because our simulation represents a special scenario where the density is entirely flat along the ridge, making the ridge a subset of critical points where ∇f = 0. Hence, the convergence enjoys a faster rate than the one we develop here for a more generic setting, because the rate of convergence for gradient estimation is faster than that for Hessian estimation using KDE. Also see [29], where both scenarios are discussed for static ridge estimation.

Figure 4: Hausdorff distance between the estimated ridge R(b p) and the population ridge R(p) across varying sample sizes n. The plot axes are visualized on a log-log scale.

8

Discussion

In this paper, we resolved an important open problem in the literature surrounding Subspace Constrained Mean Shift (SCMS) algorithms by introducing a paradigm shift from static density ridges to stable density ridges. We proved that the SCMS algorithm estimates the latter, which is defined by the Jacobian of the vector field associated with the SCMS algorithm and its corresponding eigenvalues. By formalizing the ridge-regular function class, we established the uniform R-linear convergence and topological surjectivity of both the continuous flow and the discrete algorithmic maps. Furthermore, our joint analysis of statistical and computational errors exposed a computational bottleneck in the original SCMS algorithm: coupling the step size to the smoothing bandwidth forces the iteration complexity to blow up polynomially as the sample size increases. By decoupling the step size from the bandwidth, we showed that practitioners can recover the log-density ridge in O(log n) iterations, making the generalized SCMS algorithm highly feasible for massive datasets. 20

Currently, our theoretical framework is localized, focusing on the convergence and surjectivity of trajectories initiated within a compact neighborhood of the stable ridge. An open question is characterizing the behavior of the continuous vector field ξ and the SCMS algorithms beyond this neighborhood, or more broadly, over the entire space Rd . As our numerical tests in Section 7 have demonstrated, global convergence from far outside this neighborhood is empirically viable. Naturally, one can extend the SCMS flow backwards in time to define the basins of attraction for the stable ridge. In classical Morse theory [24], the stable manifolds (basins of attraction) of local maxima for a Morse function partition the domain, covering the entire space up to a set of Lebesgue measure zero corresponding to the boundaries between adjacent basins. See, e.g., [4], for the analysis of hill-climbing (including Mean Shift) algorithms for the purpose of space partitioning (or clustering). Whether the attraction regions of stable ridges under the SCMS flow similarly partition Rd almost everywhere remains to be investigated. Characterizing the geometric and topological conditions on the density function f under which this global partitioning occurs represents an interesting direction for future research.

References [1] Arias-Castro, E., D. Mason, and B. Pelletier (2016). Errata: On the estimation of the gradient lines of a density and the consistency of the mean-shift algorithm. Journal of Machine Learning Research 17 (206), 1–4. [2] Arias-Castro, E. and W. Qiao (2023a). Moving up the cluster tree with the gradient flow. SIAM Journal on Mathematics of Data Science 5 (2), 400–421. [3] Arias-Castro, E. and W. Qiao (2023b). A unifying view of modal clustering. Information and Inference: A Journal of the IMA 12 (2), 897–920. [4] Arias-Castro, E. and W. Qiao (2025). Clustering by hill-climbing: Consistency results. The Annals of Statistics 53 (6), 2536–2562. [5] Chacón, J. E. (2015). A population background for nonparametric density-based clustering. Statistical Science 30 (4), 518–532. [6] Chen, Y., C. Genovese, and L. Wasserman (2015). Asymptotic theory for density ridges. Annals of Statistics 43 (5), 1896–1928. [7] Chen, Y.-C., C. R. Genovese, and L. Wasserman (2017). Statistical inference using the morsesmale complex. Electronic Journal of Statistics 11, 1390–1433. [8] Chen, Y.-C., S. Ho, P. E. Freeman, C. R. Genovese, and L. Wasserman (2015). Cosmic web reconstruction through density ridges: method and algorithm. Monthly Notices of the Royal Astronomical Society 454 (1), 1140–1156. [9] Cheng, Y. (1995). Mean shift, mode seeking, and clustering. IEEE Transactions on Pattern Analysis and Machine Intelligence 17 (8), 790–799. [10] Comaniciu, D. and P. Meer (2002). Mean shift: A robust approach toward feature space analysis. IEEE Transactions on Pattern Analysis and Machine Intelligence 24 (5), 603–619. [11] Eberly, D. (1996). Ridges in Image and Data Analysis. Kluwer Academic Publishers. [12] Einmahl, U. and D. M. Mason (2005). Uniform in bandwidth consistency of kernel-type function estimators. The Annals of Statistics 33 (3), 1380–1403. [13] Elsner, L. (1985). An optimal bound for the spectral variation of two matrices. Linear algebra and its applications 71, 77–80. [14] Fefferman, C., S. Ivanov, Y. Kurylev, M. Lassas, and H. Narayanan (2018). Fitting a putative manifold to noisy data. In Conference On Learning Theory, pp. 688–720. PMLR.

21

[15] Fukunaga, K. and L. Hostetler (1975). The estimation of the gradient of a density function, with applications in pattern recognition. IEEE Transactions on Information Theory 21 (1), 32–40. [16] Genovese, C. R., M. Perone Pacifico, I. Verdinelli, L. Wasserman, et al. (2014). Nonparametric ridge estimation. The Annals of Statistics 42 (4), 1511–1545. [17] Ghassabeh, Y. A., T. Linder, and G. Takahara (2013). On some convergence properties of the subspace constrained mean shift. Pattern Recognition 46 (11), 3140–3147. [18] Giné, E. and A. Guillou (2002). Rates of strong uniform consistency for multivariate kernel density estimators. Annales de l’Institut Henri Poincaré, Probabilités et Statistiques 38 (6), 907– 921. [19] Hastie, T. and W. Stuetzle (1989). Principal curves. Journal of the American Statistical Association 84 (406), 502–516. [20] Kégl, B., A. Krzyzak, T. Linder, and K. Zeger (2000). Learning and design of principal curves. IEEE transactions on pattern analysis and machine intelligence 22 (3), 281–297. [21] Kuznetsov, Y. A. Elements of applied bifurcation theory. Springer. [22] Lee, J. (2012). Introduction to Smooth Manifolds. Springer Science & Business Media. [23] Li, W. and S. Ghosal (2020). Posterior contraction and credible sets for filaments of regression functions. Electronic Journal of Statistics 14, 1707–1743. [24] Milnor, J. W. (1963). Morse theory. Number 51. Princeton university press. [25] Moews, B., J. R. Argueta Jr, and A. Gieschen (2021). Filaments of crime: Informing policing via thresholded ridge estimation. Decision Support Systems 144, 113518. [26] Ozertem, U. and D. Erdogmus (2011). Locally defined principal curves and surfaces. Journal of Machine Learning Research 12 (34), 1249–1286. [27] Perko, L. (2001). Differential Equations and Dynamical Systems (3rd ed.). Texts in Applied Mathematics. Springer-Verlag New York. [28] Qiao, W. (2021). Asymptotic confidence regions for density ridges. Bernoulli 27 (2), 946–975. [29] Qiao, W. (2025). Confidence regions for filamentary structures. Information and Inference: A Journal of the IMA 14 (4), iaaf030. [30] Qiao, W. and W. Polonik (2016). Theoretical analysis of nonparametric filament estimation. The Annals of Statistics, 1269–1297. [31] Qiao, W. and W. Polonik (2025). Algorithms for ridge estimation with convergence guarantees. Journal of Machine Learning Research 26. [32] Qiao, W. and A. Shehu (2022). Space partitioning and regression maxima seeking via a meanshift-inspired algorithm. Electronic Journal of Statistics 16 (2), 5623–5658. [33] Roweis, S. T. and L. K. Saul (2000). Nonlinear dimensionality reduction by locally linear embedding. science 290 (5500), 2323–2326. [34] Sandilya, S. and S. R. Kulkarni (2002). Principal curves with bounded turn. IEEE Transactions on Information Theory 48 (10), 2789–2793. [35] Tenenbaum, J. B., V. d. Silva, and J. C. Langford (2000). A global geometric framework for nonlinear dimensionality reduction. science 290 (5500), 2319–2323. [36] Yao, Z., J. Su, B. Li, and S.-T. Yau (2023). Manifold fitting. arXiv preprint arXiv:2304.07680 . [37] Yao, Z. and Y. Xia (2025). Manifold fitting under unbounded noise. Journal of Machine Learning Research 26 (45), 1–55. [38] Zhang, Y. and Y.-C. Chen (2023). Linear convergence of the subspace constrained mean shift algorithm: from euclidean to directional data. Information and Inference: A Journal of the IMA 12 (1), 210–311. [39] Zhang, Y. and Y.-C. Chen (2026). Mode and ridge estimation in euclidean and directional product spaces: A mean shift approach. Journal of Computational and Graphical Statistics 35 (1), 101–110.

22

A

Static vs. Stable Ridges

As discussed in the Introduction section, the local concavity of the density (λk+1 < 0) and the dynamic stability of the continuous flow (Re(µk+1 ) < 0) are distinct conditions for the static ridge and stable ridge, respectively. Because the Jacobian Jξ incorporates the rotation of the trailing eigenspace, these two conditions only coincide under highly regular curvature conditions. In general, the static ridge and stable ridge are different structures. To quantify this distinction for a 1-dimensional ridge (k = 1) in R2 , we derive the eigenvalues and eigenvectors of the Jacobian at a point x where ξ(x) = 0. By the product rule, the Jacobian of the vector field ξ(x) = Π(x)∇f (x) is Jξ (x) = Π(x)∇2 f (x) + (∇Π(x))∇f (x).

(46)

Let v1 (x) and v2 (x) be the unit eigenvectors corresponding to the leading and trailing Hessian eigenvalues (λ1 (x) ≥ λ2 (x)), so the projection matrix is Π(x) = v2 (x)v2 (x)⊤ . Below we verify that v2 is an eigenvector of Jξ (x). To show this, we multiply both sides of (46) by v2 . Since ξ(x) = 0, the first term becomes Π(x)∇2 f (x)v2 (x) = Π(x)(λ2 (x)v2 (x)) = λ2 (x)v2 (x). The second term becomes (∇v2 Π)∇f = (∇v2 v2 )(v2⊤ ∇f ) + v2 (∇v2 v2 )⊤ ∇f = v2 ⟨∇v2 v2 , ∇f ⟩. Letting ∆rot = ⟨∇v2 v2 , ∇f ⟩, we can write
Jξ (x)v2 = λ2 + ∆rot v2 .

(47)

Equation (47) proves that v2 is an eigenvector of the Jacobian. To determine its index in the sorted eigenvalues Re(µ1 ) ≥ Re(µ2 ), we multiply both sides of (46) by v1 . For the first term, Π(x)∇2 f (x)v1 = Π(x)(λ1 v1 ) = 0, because v1 is in the null space of the projection matrix Π(x). Following the same derivation for the second term, we find Jξ (x)v1 = cv2 , where c = ⟨∇v1 v2 , ∇f ⟩. Since Jξ (x)v1 = 0v1 +cv2 and Jξ (x)v2 = 0v1 +(λ2 +∆rot )v2 , the Jacobian matrix Jξ (x) expressed in the (v1 , v2 ) coordinate basis takes the lower-triangular form   0 0 . (48) c λ2 + ∆rot Hence the eigenvalues of Jξ (x) are {0, λ2 + ∆rot }, which are the diagonal entries of the triangular matrix. According to our sorting rule Re(µ1 ) ≥ Re(µ2 ), we can write µ1 = max(0, λ2 + ∆rot ) and µ2 = min(0, λ2 + ∆rot ). By the definition of the stable ridge in (3), a point x such that ξ(x) = 0 belongs to the stable ridge if and only if λ2 (x) + ∆rot (x) < 0, where the term ∆rot quantifies how the trailing eigenspace rotates as probability mass flows toward the ridge.

A.1

Examples of Ridge Divergence and Coincidence

This subsection provides three examples comparing the static ridge (Rstatic ) and the stable ridge (Rstable ) for k = 1 and d = 2. As derived in the preceding derivation, the relation between static and stable ridges is determined by the sign of λ2 + ∆rot . These examples are visualized in Figure 5. For each example, we only focus on the ridges restricted to a compact domain D ⊂ R2 .

23

Figure 5: Three examples illustrating the difference and coincidence between static and stable ridges. Green segments indicate coincidence (Rstable ∩ Rstatic ), red indicates static-only branches that are dynamically unstable (Rstatic \ Rstable ), and blue indicates dynamic attractors located in static valleys (Rstable \ Rstatic ). Starting from grids of initial points (gray dots), the SCMS continuous flows φt (black lines) converge to the limit points (black dots) on the stable ridges, but not always on the static ridges. A.1.1

Example 1: Qiao-Polonik Pitchfork Bifurcation

We first analyze the density f (u, v) = C + 38 (1 − u2 )v restricted to D = [−0.9, 0.9] × [0.1, 1.5]. This was introduced by [31] as a counterexample to show that SCMS does not always converge to the static ridge. Here, and in the following examples, a positive constant C > 0 ensures f is positive on D and is a valid probability density, and ridges remain invariant to a constant shift in the density. Let c = 3/8. The gradient and Hessian of f are     −2uv −2v −2u 2 ∇f (u, v) = c , ∇ f (u, v) = c . 1 − u2 −2u 0 √
The eigenvalues of the Hessian are λ1,2 = c −v ± v 2 + 4u2 . Both ridge definitions require ξ(u, v) = 0, which corresponds to the union of two sets: the vertical 2 axis (u = 0) and a pair of symmetric wings defined by v = ψ(u) = √ 1−u 2 . 2(1+u )

The static ridge: The condition for the static ridge is λ2 < 0. On the vertical axis (u = 0), the trailing eigenvalue is λ2 = −2cv,√which is
negative for all v > 0. On the symmetric wings, the trailing eigenvalue is λ2 = c −v − v 2 + 4u2 , which is negative, since both v and c are positive across the domain. Therefore, the static ridge claims the entire pitchfork structure, including the vertical axis extending down to the bottom of the domain D. The stable ridge: We compute ∆rot . On the vertical axis, the trailing eigenvector is v2 = (1, 0)⊤ . It follows from straightforward calculations that ∇v2 v2 = (0, 1/v)⊤ and ∇f = (0, c)⊤ , which gives ∆rot = vc . The dynamic stability condition on the axis requires µ2 (0, v) = λ2 + ∆rot = −2cv +

c 1 < 0 =⇒ v > √ . v 2

√ Consequently, the lower axis below 1/ 2 loses dynamic stability in the stable h and is 2not 2included i p (1−u )(3u −1) 2 ridge. Conversely, by calculation we have µ2 = −c 2(1 + u ) 1 + (3u2 +1)2 along the wings 24

v = ψ(u), which remains negative for all u ̸= 0. Thus, the symmetric wings are retained in the √ stable ridge. The stable ridge drops the lower axis, defining a bifurcation at v = 1/ 2. For this example, the stable ridge is a strict subset of the static ridge. A.1.2

Example 2: Static Valley as Stable Ridge

Consider D = [−0.5, 0.5] × [−0.1, 0.4] and a density defined as 1 f (u, v) = C + v + 5v 2 + (1 + 18v)u2 . 2     u(1 + 18v) 1 + 18v 18u 2 The derivatives are ∇f (u, v) = and ∇ f (u, v) = . 1 + 10v + 9u2 18u 10 Setting ξ(u, v) = 0 reveals that the vertical axis (u = 0) is the only candidate set for both the static and stable ridges. On this axis, the Hessian is diagonal with ∇2 f (0, v) = diag(1 + 18v, 10). Throughout D where v ≤ 0.4, the trailing eigenvalue is λ2 = 1 + 18v, and the leading eigenvalue is λ1 = 10. The static ridge: The axis belongs to Rstatic only where λ2 < 0. Within D, this restricts the static ridge to the sub-interval v ∈ [−0.1, −1/18). For any v > −1/18, we have λ2 > 0. Geometrically, this positive curvature means the density is locally convex in the normal direction, making that portion of the vertical axis a static valley rather than a static ridge. The stable ridge: We compute ∆rot . The trailing eigenvector is v2 = (1, 0)⊤ and  the leading  0 18 eigenvector is v1 = (0, 1)⊤ . The directional derivative of the Hessian is ∇v2 ∇2 f = . By 18 0 spectral perturbation theory, we have     18 v ⊤ (∇v2 ∇2 f )v2 0 0 v1 = = . ∇v 2 v2 = 1 2 λ2 − λ 1 (1 + 18v) − 10 1 2v−1 and Combining this with ∇f = (0, 1 + 10v)⊤ yields ∆rot = 2(1+10v) 2v−1 λ2 + ∆rot = (1 + 18v) +

36v 2 + 4v + 1 2(1 + 10v) = . 2v − 1 2v − 1

Because 36v 2 + 4v + 1 > 0 for all v, the sign of λ2 + ∆rot is determined by the denominator 2v − 1, which is negative for all v ∈ [−0.1, 0.4], ensuring µ2 = λ2 + ∆rot < 0 on the entire vertical axis (u = 0) within D. The stable ridge spans the entire vertical line v ∈ [−0.1, 0.4], even where it is a static valley for v > −1/18. For this example, the static ridge is a strict subset of the stable ridge. A.1.3

Example 3: Coincidence Between Static and Stable Ridges

We construct an example for which the two ridge definitions coincide, √ that is, Rstable = Rstatic . Consider a radially symmetric density f (u, v) = g(ρ), where ρ = u2 + v 2 , formed by adding isotropic Gaussian noise to a uniform distribution on a 1-dimensional circle of radius r. We analyze a tubular neighborhood D where the profile g(ρ) has a global maximum at ρ = r∗ and g ′′ (ρ) < 0. Under the standard polar transformation u = ρ cos θ and v = ρ sin θ, we define the orthonormal basis vectors eρ = (cos θ, sin θ)⊤ and eθ = (− sin θ, cos θ)⊤ . (Note that while eρ indicates the direction of radial extension, its Cartesian orientation depends solely on the angle θ). The gradient is radial, that is, ∇f (u, v) = g ′ (ρ)eρ . By applying ∇eρ = ρ1 eθ e⊤ θ , the Hessian is ∇2 f (u, v) = g ′′ (ρ)eρ e⊤ ρ + 25

g ′ (ρ) eθ e⊤ θ . ρ

The trailing eigenvector is v2 = eρ , with trailing eigenvalue λ2 = g ′′ (ρ). The projected gradient is ξ(u, v) = g ′ (ρ)eρ , which vanishes if and only if g ′ (ρ) = 0, occurring at ρ = r∗ . The static ridge: Because g ′′ (r∗ ) < 0, the static ridge is just the 1-dimensional circle of radius ∗ r . The stable ridge: For this example, ∆rot = 0, because   1 1 ⊤ eθ eθ eρ = eθ (e⊤ ∇v2 v2 = (∇eρ )eρ = θ eρ ) = 0. ρ ρ Consequently, µ2 ≡ λ2 = g ′′ (r∗ ), confirming that the stable and static ridges are identical in this example.

B

Proofs

B.1

Proofs for Section 3

Proof of Lemma 1. By Assumption (A1), f ∈ C 4 (Rd , R). By Assumption (A2), the condition η > 0 ensures both Π and ξ are C 2 -smooth on int(D). Consider a point x ∈ R(f ). Let V (x) ∈ Rd×(d−k) be a matrix whose columns form an orthogonal basis for the trailing eigenspace E⊥ (x). By Property 2 of Definition 1, the symmetric part of the projected Jacobian V (x)⊤ Jξ (x)V (x) is negative definite. Thus, V (x)⊤ Jξ (x)V (x) is non-singular and has full rank d − k. Let U ⊂ int(D) be a sufficiently small neighborhood around x. We construct a C 2 -smooth map f (y) = Π(y)V (x). ψ : U → Rd−k such that ξ(y) = 0 is equivalent to ψ(y) = 0 for all y ∈ U . Define W f Because Π(x)V (x) = V (x) has full column rank d − k and Π is continuous, W (y) retains full column rank on U for sufficiently small U . Applying Gram-Schmidt orthogonalization to the columns of f (y) yields an orthonormal basis for E⊥ (y), denoted by the matrix W (y) ∈ Rd×(d−k) . Because the W projection matrix Π and the Gram-Schmidt operations are C 2 -smooth, W is C 2 -smooth on U . Note that W (x) = V (x) since the columns of V (x) are already orthonormal. By construction, Π(y) = W (y)W (y)⊤ , yielding ξ(y) = W (y)ψ(y), where ψ(y) = W (y)⊤ ∇f (y).

(49)

Since W (y) has full column rank, the condition ξ(y) = 0 is equivalent to ψ(y) = 0 for y ∈ U . Furthermore, ψ is C 2 -smooth on U since both W and ∇f are C 2 -smooth. Let ψi denote the i-th component of ψ, and Wi denote the i-th column of W for i = 1, . . . , d − k. By the multivariate product rule, the Jacobian of ξ(y) = W (y)ψ(y) evaluated at x is Jξ (x) = W (x)Jψ (x) +

d−k X

ψi (x)∇Wi (x) = V (x)Jψ (x),

(50)

i=1

because ψ(x) = 0. Left-multiplying by V (x)⊤ yields V (x)⊤ Jξ (x) = Jψ (x). Since V (x)⊤ Jξ (x)V (x) is non-singular, Jψ (x) has full row rank d − k. By the Regular Level Set Theorem [22, Corollary 5.14], ψ −1 (0) forms a C 2 -smooth surface of dimension k within U . Because x ∈ R(f ) is arbitrary, R(f ) is a global C 2 -smooth k-dimensional submanifold. Fix an ϵ ∈ (0, κ]. By Property 3 of Definition 1, any point z ∈ D satisfying ∥ξ(z)∥ ≤ ϵ ≤ κ must satisfy Re(µk+1 (z)) ≤ −γ/2. Consequently, the set can be equivalently expressed as n γo Rϵ (f ) = z ∈ D : ∥ξ(z)∥ ≤ ϵ and Re(µk+1 (z)) ≤ − . (51) 2 26

Since ξ and Re(µk+1 ) are continuous functions, Rϵ (f ) is a closed subset of the compact domain D, making it compact. Property 2 dictates Rϵ (f ) ⊂ int(D). Thus, R(f ) = ξ −1 (0) ∩ Rϵ (f ) is the intersection of a closed set and a compact set, which is compact. Finally, because R(f ) ⊂ int(D), it has no boundary. ■ Proof of Lemma 2. Let x ∈ Rϵ (f ). Using the definition in (6), we have   d d 2 ⊤ ∥ξ(φt (x))∥ = 2ξ(φt (x)) ξ(φt (x)) dt dt = 2ξ(φt (x))⊤ Jξ (φt (x))ξ(φt (x)).

(52) (53)

Since ξ(φt (x)) ∈ E⊥ (φt (x)), Property 2 of Definition 1 implies d ∥ξ(φt (x))∥2 ≤ −2γ∥ξ(φt (x))∥2 . dt

(54)

By Grönwall’s inequality, (54) implies ∥ξ(φt (x))∥2 ≤ e−2γt ∥ξ(x)∥2 =⇒ ∥ξ(φt (x))∥ = ∥φ̇t (x)∥ ≤ e−γt ∥ξ(x)∥.

(55)

Since e−γt ≤ 1 for t ≥ 0, we have ∥ξ(φt (x))∥ ≤ ∥ξ(x)∥ ≤ ϵ. Furthermore, because ϵ ≤ κ, Property 3 of Definition 1 guarantees Re(µk+1 (φt (x))) ≤ −γ/2 < 0. Therefore, φt (x) ∈ Rϵ (f ) for all t ≥ 0. The total arc length L(x) of φt (x) satisfies Z ∞ Z ∞ 1 (56) L(x) = ∥φ̇t (x)∥ dt ≤ e−γt ∥ξ(x)∥ dt = ∥ξ(x)∥. γ 0 0 Since the integral on the left-hand side converges, the trajectory φt (x) satisfies the Cauchy criterion as t → ∞ in Rd . Thus, φt (x) converges to a unique limit point Φ(x) = limt→∞ φt (x). By (55), limt→∞ ∥ξ(φt (x))∥ = 0. By the continuity of ξ on D, ∥ξ(Φ(x))∥ = 0. Because ∥ξ(Φ(x))∥ = 0 < κ, Property 3 of Definition 1 implies Re(µk+1 (Φ(x))) ≤ − γ2 < 0. Thus, Φ(x) ∈ R(f ). ■ of Proof of Corollary 1. Since R(f ) ⊆ Rϵ (f ), we have d(R(f )|Rϵ (f )) = 0, and thus dH (Rϵ (f ), R(f )) = d(Rϵ (f )|R(f )). Let x ∈ Rϵ (f ). By Lemma 2, we have Φ(x) ∈ R(f ), which yields inf

∥x − z∥ ≤ ∥x − Φ(x)∥ ≤ L(x) ≤

z∈R(f )

∥ξ(x)∥ ϵ ≤ . γ γ

(57)

Taking the supremum over x ∈ Rϵ (f ) establishes the upper bound in (11). Next we consider the lower bound. Let x̃ ∈ Rϵ (f ) such that ∥ξ(x̃)∥ = ϵ. Because R(f ) is compact by Lemma 1, there exists x∗ ∈ R(f ) such that inf z∈R(f ) ∥x̃ − z∥ = ∥x̃ − x∗ ∥. Denote v = x̃ − x∗ . By (57), we have ∥v∥ ≤ ϵ/γ. We first show the Euclidean ball Bϵ/γ (x∗ ) is contained in Rκ (f ). For any y ∈ Bϵ/γ (x∗ ), define the line segment z(t) = x∗ + t(y − x∗ ) for t ∈ [0, 1]. Because x∗ ∈ R(f ) is in the interior of Rκ (f ), the segment begins in D. As long as the partial segment z([0, t]) ⊆ D, the γ κ. This implies mean value theorem ensures ∥ξ(z(t))∥ ≤ L∥z(t) − x∗ ∥ ≤ L(ϵ/γ) ≤ κ, due to ϵ ≤ L the entire segment z([0, 1]) cannot exit Rκ (f ) ⊆ D, guaranteeing Bϵ/γ (x∗ ) is contained in D.

27

Because Bϵ/γ (x∗ ) is convex, the line segment connecting x∗ and x̃ lies entirely within D. Since ξ ∈ C 2 (D) and ξ(x∗ ) = 0, a Taylor expansion yields ξ(x̃) = ξ(x∗ ) + Jξ (x∗ )v + R(v) = Jξ (x∗ )v + R(v),

(58)

where R(v) is the second-order remainder term. By the triangle inequality, we can write ∥ξ(x̃)∥ = ∥Jξ (x∗ )v + R(v)∥ ≤ ∥Jξ (x∗ )v∥ + ∥R(v)∥ ≤ L∥v∥ + M ∥v∥2 ,

(59)

where in the last inequality we have used the definition L = supx∈D ∥Jξ (x)∥ and M = 12 supx∈D ∥∇2 ξ(x)∥. Substituting ∥ξ(x̃)∥ = ϵ into (59) yields ϵ ≤ L∥v∥ + M ∥v∥2 ≤ (L + M ϵ/γ)∥v∥,

(60)

where we use ∥v∥ ≤ ϵ/γ from (57). Because ϵ ≤ Lγ M , it follows that inf

∥x̃ − z∥ = ∥v∥ ≥

z∈R(f )

1 1 ϵ≥ ϵ, L + M ϵ/γ 2L

establishing the lower bound in (11).

(61) ■

Proof of Theorem 1. Let x∗ ∈ R(f ). By definition, ξ(x∗ ) = 0, making x∗ an equilibrium point of ξ. Because R(f ) is a k-dimensional submanifold of D (Lemma 1), Jξ (x∗ ) possesses k zero eigenvalues corresponding to the tangent space. By Property 2 of Definition 1, the remaining d − k eigenvalues have real parts upper bounded by −γ < 0. The Center Manifold Theorem [27, Sec 2.7] guarantees s the existence of a C 2 -smooth, (d − k)-dimensional manifold Wloc (x∗ ) locally passing through x∗ , defined as  s Wloc (x∗ ) = y ∈ U : lim φt (y) = x∗ , (62) t→∞

s where U is a neighborhood of x∗ . By (10), we have Φ(y) = x∗ for all y ∈ Wloc (x∗ ). s Let g(x) = ∥ξ(x)∥. Because k < d, we can choose y ∈ Wloc (x∗ ) \ {x∗ } such that 0 < g(y) < ϵ. Consider the trajectory z(t) governed by the reversed flow ż(t) = −ξ(z(t)) with initial condition z(0) = y. Following the derivation in Lemma 2, while z(t) ∈ Rϵ (f ), we have

d ∥ξ(z(t))∥2 = −2ξ(z(t))⊤ Jξ (z(t))ξ(z(t)) ≥ 2γ∥ξ(z(t))∥2 > 0. dt

(63)

By Grönwall’s inequality, we obtain ∥ξ(z(t))∥2 ≥ e2γt ∥ξ(y)∥2 . Because ∥ξ(z(t))∥ grows exponentially, z(t) cannot remain in the bounded sub-level set {x ∈ D : g(x) < ϵ} for all t ≥ 0. Since Rϵ (f ) is a compact set contained in int(D) (see (51)), the continuous trajectory z(t) must intersect the boundary ∂Rϵ (f ). By the Intermediate Value Theorem, there exists a first exit time t0 > 0 such that g(z(t0 )) = ϵ. For all t ∈ [0, t0 ], we have g(z(t)) ≤ ϵ ≤ κ. By Property 3 of Definition 1, this ensures Re(µk+1 (z(t))) ≤ −γ/2 < 0, which validates the domain of the expansion bound up to t0 . Let x0 = z(t0 ). Because g(x0 ) = ϵ and Re(µk+1 (x0 )) ≤ −γ/2 < 0, it follows that x0 ∈ ∂Rϵ (f ). Because x0 is reached by the reversed flow z(t) from y, the forward flow φt (x0 ) passes through s y ∈ Wloc (x∗ ). Thus, limt→∞ φt (x0 ) = limt→∞ φt (y) = x∗ , which yields Φ(x0 ) = x∗ . Since x∗ ∈ R(f ) is arbitrary, Φ is surjective. ■ Proof of Lemma 3. As stated in the beginning of the proof of Lemma 1, ξ is C 2 -smooth on int(D). This implies that for any finite time t ≥ 0, the flow φt (x) is continuously differentiable with respect to its initial condition x ∈ Rκ (f ) ⊂ int(D). 28

Let Yx (t) = Jφt (x). Differentiating φ̇t (x) = ξ(φt (x)) with respect to x yields the variational equation Ẏx (t) = Jξ (φt (x))Yx (t), with Yx (0) = I. Let z = Φ(x) ∈ R(f ). By defining Az = Jξ (z) and Bz (t) = Jξ (φt (x)) − Jξ (z), we can write Ẏx (t) = [Az + Bz (t)] Yx (t). Applying the variation of constants formula yields Z t Yx (t) = eAz t + eAz (t−s) Bz (s)Yx (s) ds.

(64)

0

By Property 2 of Definition 1, the spectrum of Az consists of an eigenvalue 0 of algebraic multiplicity k, and d − k eigenvalues with real parts bounded above by −γ < 0. Because R(f ) is a k-dimensional manifold of equilibria, the tangent space Tz R(f ) is contained in the null space of Az . Thus, the geometric multiplicity of the zero eigenvalue equals k. Then by the real Schur decomposition, there exists an orthogonal matrix Uz = [U1z , U2z ] such that the k columns of U1z form an orthonormal basis for the null space of Az . Since Az U1z = 0, we obtain   0k×k A1z ⊤ ⊤ Az = Uz (Uz Az Uz )Uz = Uz U ⊤, 0(d−k)×k A2z z ⊤ ⊤ where A1z = U1z Az U2z and A2z = U2z Az U2z . Here A2z ∈ R(d−k)×(d−k) is a quasi-upper-triangular matrix whose 1 × 1 and 2 × 2 diagonal blocks correspond to the non-zero eigenvalues of Az , all of which have real parts bounded above by −γ. Define Ez (t) = Uz⊤ eAz t Uz , which satisfies Ėz (t) = (Uz⊤ Az Uz )Ez (t) with Ez (0) = Ik×k . Partitioning Ez (t) into 2 × 2 blocks of the same sizes as Az yields  ′    ′  Ė1z (t) Ė1z (t) 0k×k A1z E1z (t) E1z (t) Ėz (t) = = . ′ ′ 0(d−k)×k A2z E2z (t) E2z (t) Ė2z (t) Ė2z (t) ′ ′ It can be seen that E2z (t) = 0, E1z (t) = Ik×k , and E2z (t) = eA2z t . Subsequently Ė1z (t) = A1z eA2z t , that is, Z t A2z t E1z (t) = A1z eA2z s ds = A1z A−1 − I) = A1z A−1 (65) 2z (e 2z (E2z (t) − I). 0

Recall that L = supy∈D ∥Jξ (y)∥. Since Uz is orthogonal, we have max(∥A1z ∥, ∥A2z ∥) ≤ ∥Az ∥ = ∥Jξ (z)∥ ≤ L < ∞,

(66)

for all z ∈ R(f ). Define the set of matrices n o W = W ∈ R(d−k)×(d−k) : ∥W ∥ ≤ L and max Re(λj (W )) ≤ −γ . j

By continuity of the matrix 2-norm and eigenvalues, W is compact. By definition, A2z ∈ W for all z ∈ R(f ). Define ψ : W × [0, ∞) → R by ψ(W, t) = ∥eW t ∥eγt/2 = ∥e(W +(γ/2)I)t ∥. 29

For W ∈ W, the eigenvalues of W +(γ/2)I have real parts bounded above by −γ/2 < 0. This ensures limt→∞ ψ(W, t) = 0. Then since ψ is continuous on W, ψ attains a maximum on W × [0, ∞). Let Cmax = max(1, sup(W,t) ψ(W, t)) < ∞. Then for all z ∈ R(f ), ∥E2z (t)∥ = ψ(A2z , t)e−γt/2 ≤ Cmax e−γt/2 .

(67)

−1 ∥E1z (t)∥ ≤ ∥A1z ∥∥A−1 2z ∥(∥E2z (t)∥ + 1) ≤ ∥A1z ∥∥A2z ∥(Cmax + 1).

(68)

Using (65), it follows that

Because eAz t = Uz Ez (t)Uz⊤ and Uz is orthogonal,   Ik×k E1z (t) Az t ≤ 1 + ∥E1z (t)∥ + ∥E2z (t)∥. ∥e ∥ = ∥Ez (t)∥ = 0(d−k)×k E2z (t)

(69)

By the definition of W, 0 is not an eigenvalue of any W ∈ W, and hence every matrix in W is invertible. Since the map W 7→ ∥W −1 ∥ is continuous on the compact set W, we have Kγ := maxW ∈W ∥W −1 ∥ < ∞. Therefore, ∥A−1 2z ∥ ≤ Kγ for all z ∈ R(f ). It then follows from (66), (67), (68), and (69) that ∥eAz t ∥ ≤ 1 + LKγ (Cmax + 1) + Cmax =: κmax < ∞.

(70)

Recall that M = 21 supx∈D ∥∇2 ξ(x)∥ < ∞. Since limτ →∞ φτ (x) = z, we can write Z ∞ Bz (t) = Jξ (φt (x)) − Jξ (z) = − ∇2 ξ(φτ (x))ξ(φτ (x)) dτ. t −γτ

−γτ

By Lemma 2, ∥ξ(φτ (x))∥ ≤ e

∥ξ(x)∥ ≤ κe . As a result, Z ∞ 2M κ −γt e . ∥Bz (t)∥ ≤ 2M κe−γτ dτ = γ t

(71)

Plugging the bounds in (70) and (71) into (64), we have   Z t 2M κ −γs ∥Yx (t)∥ ≤ κmax + κmax e ∥Yx (s)∥ ds. γ 0 By Grönwall’s inequality, we obtain Z t    2κmax M κ −γs 2κmax M κ ∥Yx (t)∥ ≤ κmax exp e ds ≤ κmax exp =: LΦ < ∞. γ γ2 0 It follows from (67) that limt→∞ E2z (t) = 0, which, via (65), further leads to −1 lim E1z (t) = lim A1z A−1 2z (E2z (t) − I) = −A1z A2z ,

t→∞

t→∞

and Dz (∞) := lim e t→∞

Az t

= lim Uz Ez (t)Uz⊤ = Uz t→∞

 Ik×k 0

 −A1z A−1 2z U ⊤ . z 0

In (64), for the integral term, because ∥eAz (t−s) ∥ ≤ κmax for all t ≥ s and ∥Bz (s)Yx (s)∥ decays exponentially, the integrand eAz (t−s) Bz (s)Yx (s) is pointwise convergent as t → ∞ and is uniformly 30

dominated by an integrable function. By the Dominated Convergence Theorem, the limit Yx (∞) = limt→∞ Yx (t) exists in the form of Z ∞ Yx (∞) = Dz (∞) + Dz (∞)Bz (s)Yx (s) ds. 0

and satisfies ∥Yx (∞)∥ ≤ LΦ . In fact, we can establish uniform convergence of Yx (t) as t → ∞ over x ∈ Rκ (f ). Notice that Z t Z ∞ Az t Az (t−s) Yx (t) − Yx (∞) = (e − Dz (∞)) + (e − Dz (∞))Bz (s)Yx (s) ds − Dz (∞)Bz (s)Yx (s) ds. 0

t

(72) For any τ ≥ 0, by using (65), we have    0 E1z (τ ) + A1z A−1 0 2z U ⊤ = U eAz τ − Dz (∞) = Uz z z 0 E2z (τ ) 0

 A1z A−1 2z E2z (τ ) U ⊤ , z E2z (τ )

and subsequently by (66), (67) and (68), −γτ /2 ∥eAz τ − Dz (∞)∥ ≤ (1 + ∥A1z ∥∥A−1 ≤ κmax e−γτ /2 . 2z ∥)∥E2z (τ )∥ ≤ (1 + LKγ )Cmax e

where we have used (70), and by using it again we have ∥Dz (∞)∥ ≤ 2κmax . Plugging these established bounds as well as (71) and ∥Yx (s)∥ ≤ LΦ into (72) yields ∥Yx (t) − Yx (∞)∥     Z ∞ 2M κLΦ 2M κLΦ e−γs ds + 2 κmax e−γs ds κmax e−γ(t−s)/2 γ γ t 0     Z t 2M κLΦ 2M κLΦ = κmax e−γt/2 + κmax e−γt/2 e−γt e−γs/2 ds + 2κmax 2 γ γ 0   8M κLΦ −γt/2 −γt/2 ≤ κmax e + κmax e . γ2 ≤ κmax e−γt/2 +

Z t

This confirms that Yx (t) converges uniformly to Yx (∞) over Rκ (f ). Because φt (x) converges pointwise to Φ(x) and Jφt (x) = Yx (t) converges uniformly to Yx (∞), Φ(x) is continuously differentiable with JΦ (x) = Yx (∞). Thus, supx∈Rκ (f ) ∥JΦ (x)∥ ≤ LΦ . ■

B.2

Proofs for Section 4

Proof of Lemma 4. We proceed by induction on m to prove part (a). The base case m = 0 holds by assumption. Assume xm ∈ Rϵ (f ). Because ϵ ≤ κ, we have xm ∈ Rκ (f ). By Property 2 of Definition 1, we have sup

u⊤ Jξ (xm )u ≤ −γ < 0.

(73)

u∈E⊥ (xm ),∥u∥=1

By definition, we have ∥xm+1 − xm ∥ = ∥αξ(xm )∥ ≤ αϵ < d∂ , which implies that the line segment connecting xm and xm+1 is contained in int(D). This validates the application of a Taylor expansion to ξ(xm+1 ) = ξ(xm + αξ(xm )) around xm , from which we have ∥ξ(xm+1 )∥ ≤ ∥[I + αJξ (xm )]ξ(xm )∥ + M α2 ∥ξ(xm )∥2 . 31

(74)

Since ξ(xm ) = Π(xm )∇f (xm ) ∈ E⊥ (xm ), (73) implies ∥[I + αJξ (xm )]ξ(xm )∥2 = ∥ξ(xm )∥2 + 2αξ(xm )⊤ Jξ (xm )ξ(xm ) + α2 ∥Jξ (xm )ξ(xm )∥2 ≤ (1 − 2αγ + α2 L2 )∥ξ(xm )∥2

(75)

2

< (1 − αγ)∥ξ(xm )∥ , 2 where the last inequality follows from α < Lγ2 . Because (1 − αγ) < (1 − αγ 2 ) , we have  αγ  ∥ξ(xm )∥. ∥[I + αJξ (xm )]ξ(xm )∥ ≤ 1 − 2

Substituting (76) into (74) yields     αγ αγ ∥ξ(xm+1 )∥ ≤ 1 − + M α2 ∥ξ(xm )∥ ∥ξ(xm )∥ ≤ 1 − + M α2 ϵ ∥ξ(xm )∥, 2 2

(76)

(77)

αγ γ 2 where we use the inductive hypothesis ∥ξ(xm )∥ ≤ ϵ. Because α < 4M ϵ , it follows that M α ϵ < 4 . Substituting this bound into (77) yields   αγ  αγ  αγ + ∥ξ(xm )∥ = 1 − ∥ξ(xm )∥ = ρ∥ξ(xm )∥. (78) ∥ξ(xm+1 )∥ < 1 − 2 4 4

Since ρ < 1, we have ∥ξ(xm+1 )∥ ≤ ρϵ < ϵ ≤ κ. By Property 3 of Definition 1, Re(µk+1 (xm+1 )) ≤ −γ/2 < 0. Thus, xm+1 ∈ Rϵ (f ), completing the induction for part (a). Notice that (78) establishes part (b). For part (c), by definition and (78), we have ∥xm+1 − xm ∥ = ∥αξ(xm )∥ ≤ αρm ∥ξ(x0 )∥, which leads to ∞ ∞ X X α Lα (x0 ) = ∥xm+1 − xm ∥ ≤ α∥ξ(x0 )∥ ρm = ∥ξ(x0 )∥. (79) 1 − ρ m=0 m=0 Substituting ρ = 1 − αγ 4 gives Lα (x0 ) ≤

α 4
∥ξ(x0 )∥. αγ ∥ξ(x0 )∥ = γ 1− 1− 4

(80)

This establishes part (c) and concludes the proof.

■ P∞

Proof of Theorem 2. By Lemma 4, the total discrete arc length Lα (x0 ) = m=0 ∥xm+1 − xm ∥ is finite, which implies {xm } forms a Cauchy sequence in Rd . Thus, the sequence {xm } converges to a unique limit point x∞ . Because xm ∈ Rϵ (f ) for all m, and Rϵ (f ) is a closed set, it must contain its limit points, yielding x∞ ∈ Rϵ (f ). By the continuity of ξ and Lemma 4, we have ∥ξ(x∞ )∥ = limm→∞ ∥ξ(xm )∥ = 0. Property 3 of Definition 1 guarantees that Re(µk+1 (x∞ )) ≤ −γ/2 < 0, which establishes x∞ ∈ R(f ). Using the triangle inequality and Lemma 4, we can write ∥xm − x∞ ∥ ≤

∞ X j=m

∥xj+1 − xj ∥ =

∞ X

α∥ξ(xj )∥ ≤ α∥ξ(x0 )∥

j=m

∞ X

ρj = α∥ξ(x0 )∥

j=m

ρm . 1−ρ

(81)

By substituting ∥ξ(x0 )∥ ≤ ϵ and utilizing the relation 1 − ρ = αγ 4 , we obtain ∥xm − x∞ ∥ ≤ α∥ξ(x0 )∥

ρm 4∥ξ(x0 )∥ m = ρ ≤ Cρm . αγ/4 γ 32

(82)

■ Proof of Proposition 1. Given f ∈ C 4 , the vector field ξ(x) = Π(x)∇f (x) is C 2 -smooth, ensuring the one-step iteration map Gα (x) is continuous. By Lemma 4, Gα (x) ∈ Rϵ (f ) for all x ∈ Rϵ (f ). Thus, the m-fold composition Gm α (x) is a continuous function on Rϵ (f ) for all m ≥ 0. By Theorem 2, for all x ∈ Rϵ (f ) and m ≥ 0, m ∥Gm α (x) − Φα (x)∥ ≤ Cρ .

Since the bound Cρm is independent of x and limm→∞ Cρm = 0, the sequence Gm α (x) converges uniformly to Φα (x) on Rϵ (f ). By the Uniform Limit Theorem, the uniform limit of continuous functions is continuous. Thus, Φα (x) is continuous on Rϵ (f ). ■ Proof of Proposition 2. Recall xm = Gm α (x0 ). By Theorem 2, limm→∞ xm = Φα (x0 ) ∈ R(f ). Because Φ is continuous (see Lemma 3) and acts as the identity operator on R(f ), we have limm→∞ Φ(xm ) = Φ(Φα (x0 )) = Φα (x0 ). Thus, we can write ∥Φα (x0 ) − Φ(x0 )∥ =

∞ X

[Φ(xm+1 ) − Φ(xm )] ≤

m=0

∞ X

∥Φ(xm+1 ) − Φ(xm )∥.

(83)

m=0

For the continuous flow φt defined in (6), we have Φ(xm ) = Φ(φα (xm )). By definition, xm+1 = xm +αξ(xm ), which yields ∥xm+1 −xm ∥ = α∥ξ(xm )∥ R α≤ αϵ. By Lemma 2, ∥ξ(φt (xm ))∥ ≤ ∥ξ(xm )∥ ≤ ϵ for all t ≥ 0, which leads to ∥φα (xm ) − xm ∥ ≤ 0 ∥ξ(φt (xm ))∥ dt ≤ αϵ. Thus, both xm+1 and φα (xm ) belong to the closed Euclidean ball Bαϵ (xm ). For y ∈ Bαϵ (xm ), define the line segment z(s) = xm + s(y − xm ) for s ∈ [0, 1]. Assume for contradiction that z(s) ∈ / Rκ (f ) for some s. Since z(0) = xm ∈ Rϵ (f ) and z(s) is continuous, there exists a first exit time s∗ ≤ 1 such that ∥ξ(z(s∗ ))∥ = κ. The segment z([0, s∗ ]) is contained in Rκ (f ) ⊂ int(D). Applying the Mean Value Theorem along this segment yields ∥ξ(z(s∗ )) − ξ(xm )∥ ≤ L∥z(s∗ ) − xm ∥ ≤ L∥y − xm ∥ ≤ Lαϵ.

(84)

By the triangle inequality, ∥ξ(z(s∗ ))∥ ≤ ∥ξ(xm )∥ + Lαϵ ≤ ϵ + Lαϵ. Because α < κ−ϵ Lϵ , we obtain ∥ξ(z(s∗ ))∥ < κ, which contradicts ∥ξ(z(s∗ ))∥ = κ. Thus, no such s∗ exists, yielding Bαϵ (xm ) ⊂ Rκ (f ). Because Bαϵ (xm ) is convex, the line segment connecting xm+1 and φα (xm ) is contained in Bαϵ (xm ) ⊂ Rκ (f ). By Lemma 3, ∥JΦ (x)∥ ≤ LΦ for all x ∈ Rκ (f ). Applying the Mean Value Theorem along this segment, we have ∥Φ(xm+1 ) − Φ(xm )∥ = ∥Φ(xm+1 ) − Φ(φα (xm ))∥ ≤ LΦ ∥xm+1 − φα (xm )∥.

(85)

Since φ̇t = ξ(φt ), φ̈t = Jξ (φt )ξ(φt ) and ∥ξ(φt (xm ))∥ ≤ ∥ξ(xm )∥ by Lemma 2, we obtain ∥xm+1 − φα (xm )∥ = ∥φα (xm ) − (xm + αξ(xm ))∥ 1 1 1 ≤ α2 sup ∥φ̈t ∥ ≤ α2 L sup ∥ξ(φt (xm ))∥ ≤ Lα2 ∥ξ(xm )∥. 2 t∈[0,α] 2 2 t∈[0,α]

(86)

Substituting (86) into (85) and applying the result to (83) yields ∥Φα (x0 ) − Φ(x0 )∥ ≤

∞ X 1 m=0

2

LΦ Lα2 ∥ξ(xm )∥ =

∞ X 1 1 LΦ Lα α∥ξ(xm )∥ = LΦ LαLα (x0 ). 2 2 m=0

33

(87)

By Lemma 4, it follows that 1 ∥Φα (x0 ) − Φ(x0 )∥ ≤ LΦ Lα 2



   4 2LΦ L ∥ξ(x0 )∥ = ∥ξ(x0 )∥ α, γ γ

which concludes the proof.

(88) ■

Proof of Theorem 3. For Gα defined in (14), its Jacobian is JGα (x) = I + αJξ (x), x ∈ D. Because α < 1/L, we have ∥αJξ (x)∥ ≤ αL < 1 for all x ∈ D. Thus, JGα (x) is non-singular on D. By the Inverse Function Theorem, Gα is a C 1 -smooth local diffeomorphism on Rϵ (f ). Let x∗ ∈ R(f ). Because ξ(x∗ ) = 0, we have Gα (x∗ ) = x∗ . The Jacobian of Gα at x∗ is JGα (x∗ ) = I + αJξ (x∗ ). Since Jξ (x∗ ) has k zero eigenvalues corresponding to the tangent space, JGα (x∗ ) has k eigenvalues equal to 1. Let µ be one of the d−k non-zero eigenvalues of Jξ (x∗ ). The corresponding eigenvalue of JGα (x∗ ) is 1 + αµ. Because Re(µ) ≤ −γ and |µ| ≤ L, the squared modulus satisfies  αγ 2 , (89) |1 + αµ|2 ≤ 1 − 2αγ + α2 L2 < 1 − αγ < 1 − 4 where we use α2 L2 < αγ due to α < Lγ2 . In other words, JGα (x∗ ) has k eigenvalues equal to 1 and the remaining d − k eigenvalues have moduli bounded from above by ρ = 1 − αγ 4 < 1. Then the Center Manifold Theorem for discrete-time systems [21, Section 5.1.2] guarantees the existence of a C 1 -smooth, (d − k)-dimensional manifold WGs α (x∗ ) passing through x∗ . For a local neighborhood U of x∗ , this manifold is defined as n o WGs α (x∗ ) = y ∈ U : lim Gm (90) α (y) = x∗ . m→∞

By definition of Φα , we have Φα (y) = x∗ for all y ∈ WGs α (x∗ ). Let g(x) = ∥ξ(x)∥. Because k < d, we can choose y ∈ WGs α (x∗ ) \ {x∗ } such that 0 < g(y) < ϵ. Since Gα is a local diffeomorphism on Rϵ (f ), its local inverse G−1 α exists. Consider the sequence of (y ) with y = y. backward iterates defined by y−(m+1) = G−1 −m 0 α Following the derivation in Lemma 4, while y−(m+1) ∈ Rϵ (f ), we have g(y−m ) = g(Gα (y−(m+1) )) ≤ ρg(y−(m+1) ),

(91)

or equivalently, g(y−(m+1) ) ≥ ρ−1 g(y−m ). By induction, we obtain g(y−m ) ≥ ρ−m g(y). Because ρ < 1, g(y−m ) grows exponentially. Thus, the sequence {y−m } cannot remain in the sub-level set {x ∈ D : g(x) ≤ ϵ} for all m ≥ 0. There exists a finite integer M > 0 such that xout = y−M ∈ WGs α (x∗ ) satisfies g(xout ) > ϵ. Since WGs α (x∗ ) is diffeomorphic to an open ball in Rd−k , it is path-connected. Let σ : [0, 1] → s WGα (x∗ ) be a continuous path satisfying σ(0) = x∗ and σ(1) = xout . Because g(σ(s)) is continuous with g(σ(0)) = 0 and g(σ(1)) > ϵ, the Intermediate Value Theorem guarantees the existence of s0 ∈ (0, 1) such that x0 = σ(s0 ) ∈ WGs α (x∗ ) satisfies g(x0 ) = ϵ. Because g(x0 ) = ϵ ≤ κ, Property 3 of Definition 1 implies Re(µk+1 (x0 )) ≤ −γ/2 < 0, yielding x0 ∈ ∂Rϵ (f ). Because x0 ∈ WGs α (x∗ ), it follows that Φα (x0 ) = x∗ . Because x∗ ∈ R(f ) is arbitrary, Φα is surjective. ■

B.3

Proofs for Section 5

Proof of Lemma 6. Assume the event En (c) holds. Define δ = η/4. Since the density derivative convergence rates satisfy r1 (n) = o(r2 (n)) and r2 (n) = o(r3 (n)), there exists an integer N0 such that for all n ≥ N0 , we have r1 (n) ≤ r2 (n) ≤ r3 (n) and dC2 r2 (n) ≤ δ. 34

Part (a): By Lemma 5, the event En (c) guarantees ∥∇2 fb(x) − ∇2 f (x)∥ ≤ d max |Da fb(x) − Da f (x)| ≤ dC2 r2 (n), |a|=2

which further implies via Weyl’s inequality that bj (x) − λj (x)| ≤ ∥∇2 fb(x) − ∇2 f (x)∥ ≤ dC2 r2 (n). max |λ j

Part (b): Because supx∈D [λk (x) − λk+1 (x)] ≥ η > 0, the Davis-Kahan sin Θ theorem yields b ∥Π(x) − Π(x)∥ ≤

2dC2 2 2b ∥∇ f (x) − ∇2 f (x)∥ ≤ r2 (n), η η

for all x ∈ D.

b it follows from the triangle inequality that Part (c): Using the definition of ξ and ξ, b − ξ(x)∥ ≤ ∥Π(x)(∇ b b ∥ξ(x) fb(x) − ∇f (x))∥ + ∥(Π(x) − Π(x))∇f (x)∥. b Because ∥Π(x)∥ = 1 and Assumption A1 guarantees the existence of a constant b1 > 0 such that ∥∇f (x)∥ ≤ b1 for all x ∈ D, we obtain   √ b − ξ(x)∥ ≤ dC1 r1 (n) + 2dC2 r2 (n) b1 . ∥ξ(x) η √ We define the constant Cξ = dC1 + 2dCη2 b1 . Because n ≥ N0 ensures r1 (n) ≤ r2 (n), we have b − ξ(x)∥ ≤ Cξ r2 (n). ∥ξ(x) Part (d): The Jacobian of ξ(x) = Π(x)∇f (x) is given by Jξ (x) = Π(x)∇2 f (x) + T (∇Π(x), ∇f (x)), where the tensor operator T produces a d × d matrix with entries [T (∇Π, ∇f )]i,k = Applying the triangle inequality, we get

∂Πi,j ∂f j=1 ∂xk ∂xj .

Pd

b 2 fb − Π∇2 f ∥ + ∥T (∇Π, b ∇fb) − T (∇Π, ∇f )∥ ∥Jξb(x) − Jξ (x)∥ ≤ ∥Π∇ b 2 fb − Π∇2 f ∥ + ∥T (∇Π b − ∇Π, ∇fb)∥ + ∥T (∇Π, ∇fb − ∇f )∥. ≤ ∥Π∇

(92)

Because D is compact and f ∈ C 4 , ∥∇Π∥ and ∥∇f ∥ are uniformly bounded on D. Because ∇fb converges uniformly to ∇f by Lemma 5, it is uniformly bounded for all n ≥ N0 . Let Sd denote the space of d×d symmetric matrices, and let P : Sd → Rd×d denote the projection mapping onto the subspace spanned by the top k eigenvectors, such that Π(x) = P(∇2 f (x)) and b Π(x) = P(∇2 fb(x)). Because λk − λk+1 is uniformly bounded below by η, the image set ∇2 f (D) is a compact subset of Sd . Define the compact δ-tube Uδ = {A ∈ Sd : inf x∈D ∥A − ∇2 f (x)∥ ≤ δ}. Because δ = η/4, Weyl’s inequality guarantees every matrix in Uδ has a gap between the kth and (k + 1)th eigenvalues strictly bounded below by η/2 > 0. On this domain, the first and second differentials of P, denoted DP and D2 P, exist and are uniformly bounded in operator norm by finite constants cη and Lη on Uδ . Applying the chain rule, the partial derivative of Π with respect to the i-th coordinate xi is ∂xi Π(x) = DP(∇2 f (x))[∂xi ∇2 f (x)]. Using the triangle inequality, we have b ∥∂xi Π(x) − ∂xi Π(x)∥ ≤ ∥DP(∇2 fb(x)) − DP(∇2 f (x))∥∥∂xi ∇2 fb(x)∥ + ∥DP(∇2 f (x))∥∥∂x ∇2 fb(x) − ∂x ∇2 f (x)∥. i

35

i

(93)

Since n ≥ N0 is large enough that ∥∇2 fb(x)−∇2 f (x)∥ ≤ δ for all x ∈ D, the straight line segment connecting ∇2 f (x) and ∇2 fb(x) lies within a closed, convex δ-ball in Sd centered at ∇2 f (x), which is contained in Uδ . Applying the Mean Value Theorem along this convex segment yields ∥DP(∇2 fb(x)) − DP(∇2 f (x))∥ ≤ sup ∥D2 P(A)∥∥∇2 fb(x) − ∇2 f (x)∥ ≤ Lη (dC2 r2 (n)). A∈Uδ

Bounding the remaining terms in (93) we obtain
b ∥∂xi Π(x) − ∂xi Π(x)∥ ≤ Lη (dC2 r2 (n)) ∥∂xi ∇2 f (x)∥ + d2 C3 r3 (n) + cη (d2 C3 r3 (n)). Given f ∈ C 4 , ∥∂xi ∇2 f (x)∥ is uniformly bounded on D. Therefore, there exists a constant CΠ′ > 0 b such that supx∈D ∥∇Π(x) − ∇Π(x)∥ ≤ CΠ′ r3 (n) by recalling that n ≥ N0 ensures r2 (n) ≤ r3 (n). Substituting these bounds into (92) establishes the existence of a constant CJ > 0 satisfying sup ∥Jξb(x) − Jξ (x)∥ ≤ CJ r3 (n).

x∈D

■ Proof Proofs of Corollary 2. Let event En (c) hold. Recall L = supz∈D ∥Jξ (z)∥. Define ϵΠ (n) = 2dC2 2 1−1/d . Because rℓ (n) → 0 for η r2 (n), ϵQ (n) = 2γϵΠ (n) + 2LϵΠ (n) + LϵΠ (n) , and Cµ = (2L + 1) ℓ = 0, . . . , 4, and ϵξ (n), ϵJ (n), ϵQ (n) → 0, we can choose N1 ≥ N0 such that for all n ≥ N1 ,  d ! γ η κ ϵJ (n) ≤ min 1, , ϵξ (n) < , 2dC2 r2 (n) ≤ , 5Cµ 2 2 ϵQ (n) + ϵJ (n) ≤

γ , 2

max Cℓ rℓ (n) ≤ 1.

0≤ℓ≤4

Part (a): By Assumption (B1), we have K ∈ C 4 (Rd , R), which implies fb ∈ C 4 (Rd , R). Under Assumption (A1), max|a|=ℓ supx∈D |Da f (x)| ≤ bℓ for some constants bℓ > 0, ℓ = 0, 1, 2, 3, 4. By Lemma 5 and the triangle inequality, we have max|a|=ℓ supx∈D |Da fb(x)| ≤ bℓ +Cℓ rℓ (n) ≤ bℓ +1 =: Bℓ , because n ≥ N1 ensures Cℓ rℓ (n) ≤ 1. Part (b): By Lemma 6, part (a), for all x ∈ D, bk (x) − λ bk+1 (x) ≥ η − 2dC2 r2 (n). λ bk (x) − λ bk+1 (x) ≥ η/2. For n ≥ N1 , 2dC2 r2 (n) ≤ η/2, leading to λ b b Part (c): Let x ∈ Rκ/2 (f ). By definition, ∥ξ(x)∥ ≤ κ/2. By the triangle inequality and Lemma 6, part (c), κ κ b ∥ξ(x)∥ ≤ ∥ξ(x)∥ + ϵξ (n) < + = κ, 2 2 which implies Re(µk+1 (x)) ≤ −γ/2, by Property 3 of Definition 1. Because ∥ξ(x)∥ ≤ κ and Re(µk+1 (x)) ≤ −γ/2, we obtain x ∈ Rκ (f ). b⊥ (x) with ∥u∥ = 1, which guarantees Π(x)u b Let u ∈ E = u. Decompose u = v + w for v = Π(x)u ∈ E⊥ (x) and w = (I − Π(x))u ∈ E∥ (x), such that 1 = ∥u∥ ≤ ∥v∥ + ∥w∥. By Lemma 6, part (b), b b ∥w∥ = ∥u − Π(x)u∥ = ∥Π(x)u − Π(x)u∥ ≤ ∥Π(x) − Π(x)∥∥u∥ ≤ ϵΠ (n).

36

Thus, 1 − ϵΠ (n) ≤ 1 − ∥w∥ ≤ ∥v∥ ≤ 1. It is clear that u⊤ Jξ (x)u = v ⊤ Jξ (x)v + 2v ⊤ Jξ (x)w + w⊤ Jξ (x)w.

(94)

Because v ∈ E⊥ (x) and x ∈ Rκ (f ), it follows from Property 2 of Definition 1 that v ⊤ Jξ (x)v ≤ −γ∥v∥2 ≤ −γ(1 − ϵΠ (n))2 ≤ −γ + 2γϵΠ (n). By the Cauchy-Schwarz inequality, 2|v ⊤ Jξ (x)w| ≤ 2∥v∥∥Jξ (x)∥∥w∥ ≤ 2LϵΠ (n), |w⊤ Jξ (x)w| ≤ ∥w∥∥Jξ (x)∥∥w∥ ≤ LϵΠ (n)2 . Substituting these components into (94) yields u⊤ Jξ (x)u ≤ −γ + 2γϵΠ (n) + 2LϵΠ (n) + LϵΠ (n)2 = −γ + ϵQ (n). Then by Lemma 6, part (d), u⊤ Jξb(x)u = u⊤ Jξ (x)u + u⊤ (Jξb(x) − Jξ (x))u ≤ −γ + ϵQ (n) + ϵJ (n). b⊥ (x) we Since n ≥ N1 is large enough that ϵQ (n) + ϵJ (n) ≤ γ/2, taking the supremum over u ∈ E get the desired result in part (c). b Part (d): Suppose x ∈ D satisfies ∥ξ(x)∥ < κ/2. By Lemma 6, part (c), b b − ξ(x)∥ < ∥ξ(x)∥ ≤ ∥ξ(x)∥ + ∥ξ(x)

κ κ + = κ. 2 2

By Property 3 of Definition 1, ∥ξ(x)∥ < κ implies Re(µk+1 (x)) ≤ −γ/2. By definition, the sequences {µj (x)}dj=1 and {b µj (x)}dj=1 are sorted in descending order of their real parts. By Elsner’s Theorem [13], there exists a bijection π on the indices {1, . . . , d} such that max |b µπ(j) (x) − µj (x)| ≤ ∥Jξb(x)∥ + ∥Jξ (x)∥ j

1−1/d

∥Jξb(x) − Jξ (x)∥1/d .

By Lemma 6, part (d), we have ∥Jξb(x) − Jξ (x)∥ ≤ ϵJ (n). Because n ≥ N1 ensures ϵJ (n) ≤ 1, we have ∥Jξb(x)∥ ≤ L + 1. Thus, max |b µπ(j) (x) − µj (x)| ≤ (2L + 1)1−1/d ϵJ (n)1/d = Cµ ϵJ (n)1/d ≤ γ/5 < γ/4. j

For j ≤ k, Re(µj (x)) = 0 implies Re(b µπ(j) (x)) ≥ −γ/5. For j > k, Re(µj (x)) ≤ −γ/2 implies Re(b µπ(j) (x)) ≤ −3γ/10. Hence, {π(j) : j = 1, . . . , k} = {1, . . . , k} and {π(j) : j = k + 1, . . . , d} = {k + 1, . . . , d}. Consequently, there exists some j ′ ≥ k + 1 such that |b µk+1 (x) − µj ′ (x)| ≤ γ/5, and Re(b µk+1 (x)) ≤ Re(µj ′ (x)) +

γ γ γ γ 3γ γ ≤ Re(µk+1 (x)) + ≤ − + = − <− . 5 5 2 5 10 4

b We have shown that ∥ξ(x)∥ < κ/2 implies Re(b µk+1 (x)) < −γ/4 for x ∈ D. Equivalently, this means b Re(b µk+1 (x)) ≥ −γ/4 implies ∥ξ(x)∥ ≥ κ/2, concluding the proof of part (d). ■

37

Proof of Theorem 4. Let event En (c) hold. By Assumption (A1), f ∈ F(η, γ, κ, D). Recall L = supz∈D ∥Jξ (z)∥ and M = 12 supz∈D ∥∇2 ξ(z)∥. By Lemma 5, we can choose N2 ≥ N1 sufb := supz∈D ∥J b(z)∥ satisfies L/2 ≤ L b ≤ 2L, and ficiently large such that for all n ≥ N2 , L ξ b c := 1 supz∈D ∥∇2 ξ(z)∥ M ≤ 2M . 2 Furthermore, let N2 be large enough that   κ Lγ γκ ϵξ (n) ≤ min , , . (95) 2 8M 8L
Because n ≥ N2 ≥ N1 , Corollary 2 establishes fb ∈ F η2 , γ2 , κ2 , D . b x) = 0, which implies ∥ξ(b b x)∥ ≤ ϵξ (n). Because Let x b ∈ R(fb). By definition, ξ(b x)∥ = ∥ξ(b x) − ξ(b ϵξ (n) ≤ κ/2 < κ, Property 3 of Definition 1 implies Re(µk+1 (b x)) ≤ −γ/2. Thus, x b ∈ Rϵξ (n) (f ), establishing R(fb) ⊆ Rϵξ (n) (f ). Using (95), we can apply Corollary 1 to obtain d(R(fb)|R(f )) ≤ d(Rϵξ (n) (f )|R(f )) ≤

ϵξ (n) . γ

(96)

b Let x ∈ R(f ). By definition, we have ξ(x) = 0, and hence ∥ξ(x)∥ ≤ ϵξ (n). Because ϵξ (n) ≤ κ/2, Corollary 2 implies Re(b µk+1 (x)) ≤ −γ/4. Thus, x ∈ Rϵξ (n) (fb), establishing R(f ) ⊆ Rϵξ (n) (fb). By our choice of N2 , we have !     b (γ/2) κ κ (L/2)(γ/2) γ/2 κ κ Lγ γκ κ L(γ/2) , , ≥ min , , = min , , ≥ ϵξ (n). (97) min b 2 c 2 2 2M 2L 2 2 8M 8L L M This allows to apply Corollary 1 to fb, resulting in 2ϵξ (n) ϵξ (n) = . d(R(f )|R(fb)) ≤ d(Rϵξ (n) (fb)|R(fb)) ≤ γ/2 γ Combining the bounds from (96) and (98) yields   2ϵξ (n) ϵξ (n) 2ϵξ (n) b dH (R(f ), R(f )) ≤ max , = . γ γ γ Substituting ϵξ (n) = O(r2 (n)) establishes the final asymptotic rate.

(98)

(99) ■

Proof of Theorem 5. Assume the event En (c) holds for n ≥ N2 . By Theorem 4, N2 ensures b ≤ 2L, M c ≤ 2M , and ϵξ (n) ≤ κ/2. L Let x ∈ Rϵ (fb). By the triangle inequality and Lemma 6, b ∥ξ(x)∥ ≤ ∥ξ(x)∥ + ϵξ (n) ≤ ϵ + ϵξ (n).

(100)

Since ϵ ≤ κ/2 and ϵξ (n) ≤ κ/2, we obtain ∥ξ(x)∥ ≤ κ. This establishes Rϵ (fb) ⊆ Rκ (f ). Consequently, inf ∥x − y∥ ≥ inf ∥x − y∥ = d∂ . (101) x∈Rκ (f ),y∈∂D

x∈Rϵ (fb),y∈∂D

By Corollary 2, fb satisfies Assumption (A1) and fb ∈ F η2 , γ2 , κ2 , D . On the other hand, we can show the step size α satisfies the conditions in Theorem 2 when f is replaced by fb, that is, α<

γ γ/2 ≤ , 2 b2 8L 2L

α<

γ γ/2 ≤ , cϵ 16M ϵ 4M 38

and

α<

d∂ . ϵ

(102)

Therefore, by invoking Theorem 2, we can show the convergence of {b xm }∞ m=0 to a unique limit b point x b∞ ∈ R(f ) with the uniform bound e ρem , ∥b xm − x b∞ ∥ ≤ C

(103)

e = 4ϵ = 8ϵ and ρe = 1 − α(γ/2) = 1 − αγ . where C γ/2 γ 4 8

■

Proof of Theorem 6. Assume the event En (c) holds for n ≥ N2 . By Theorem 4, N2 ensures
b ≤ 2L, M c ≤ 2M , and ϵξ (n) ≤ κ/2. By Corollary 2, we have fb ∈ F η , γ , κ , D . L 2 2 2 As established in Theorem 5, the step size α satisfies the bounds in (102) for the empirical b b α (x) = x + αξ(x) b and M c. Furthermore, the step size satisfies α < 1 ≤ 1 , ensuring G parameters L b 2L L is a local diffeomorphism on D. b we conclude that that the limit map Φ b α (x) is surjective Applying Theorem 3 to the vector field ξ, b b from ∂Rϵ (f ) onto R(f ). ■ e = 8ϵ and ρe = 1 − αγ . Proof of Corollary 3. Let C γ 8 b α (b bm b α (b e ρem . Let x b0 ∈ ∂Rϵ (fb). By Theorem 5, we have Φ x0 ) ∈ R(fb) and ∥G x0 ) − Φ x0 )∥ ≤ C α (b Hence, bm b α (b e ρem . d(Xbm |R(fb)) ≤ sup ∥G x0 ) − Φ x0 )∥ ≤ C (104) α (b x b0 ∈∂Rϵ (fb)

b α is surjective from ∂Rϵ (fb) onto R(fb), that is, for every x∗ ∈ R(fb), Following from Theorem 6, Φ e ρem . b α (b bm x0 ) − x∗ ∥ ≤ C there exists x b0 ≡ x b0 (x∗ ) ∈ ∂Rϵ (fb) such that Φ x0 ) = x∗ . By Theorem 5, ∥G α (b Hence bm e ρem . d(R(fb)|Xbm ) ≤ sup ∥G x0 ) − x∗ ∥ ≤ C (105) α (b x∗ ∈R(fb)

e ρem . Combining (104) and (105) establishes dH (Xbm , R(fb)) ≤ C

■

Proof of Theorem 7. Let n ≥ N2 . By Lemma 5, the event En (c) occurs with probability at least 1 − n−c . Assume En (c) holds. By the triangle inequality for the Hausdorff metric, dH (Xbm , R(f )) ≤ dH (Xbm , R(fb)) + dH (R(fb), R(f )).

(106)

αγ m By Corollary 3, we have dH (Xbm , R(fb)) ≤ 8ϵ . By Theorem 4, we have dH (R(fb), R(f )) ≤ γ 1− 8 2ϵξ (n) . Substituting these bounds into (106) we obtain the final inequality. γ

B.4

■

Proofs for Section 6

Proof of Lemma 7. Let Llog = supx∈D ∥Jξlog (x)∥ and Mlog = 21 supx∈D ∥∇2 ξ log (x)∥. Under Assumptions (A1) and (A2′ ), we have p ∈ C 4 , Llog < ∞ and Mlog < ∞. Also, similar to d∂, , we define d∂,log = inf x∈Rκlog (p),y∈∂D ∥x − y∥ > 0. Let N3 be sufficiently large such that for all n ≥ N3 , conditionally on En (c), it satisfies that ϵlog ξ (n) ≤ κlog /2, ∆n ≤ 1/2, and 2

2CK h < min

γlog 1 2d∂,log γlog , , , 2 8Llog 8Mlog κlog 2Llog κlog

! ,

and ∥∇αn (x)∥ ≤ 2κ1log uniformly on D. Assume the event En (c) holds for n ≥ N3 . 39

(107)

(x) (x) The adaptive step size is αn (x) = CK h2 fqb(x) . Since ∆n = supx∈D fqb(x) − 1 ≤ 1/2, the step size satisfies the uniform bounds b

b

0<

CK h2 ≤ CK h2 (1 − ∆n ) ≤ αn (x) ≤ CK h2 (1 + ∆n ) ≤ 2CK h2 . 2

(108)

By (107), the upper bound 2CK h2 satisfies the step-size conditions of Theorem 6. b log ≡ G b log at x is The Jacobian of G αn JGblog (x) = I + αn (x)Jξblog (x) + ξblog (x)∇αn (x)⊤ .

(109)

For any x ∈ Rϵ (b p), ∥αn (x)Jξblog (x)∥ ≤ 2CK h2 Llog ≤ 1/2. In addition, since ∥ξblog (x)∥ ≤ ϵ ≤ κlog /2, we have ∥ξblog (x)∇αn (x)⊤ ∥ ≤ (κlog /2)∥∇αn (x)∥ ≤ 1/4, by the definition of N3 . These b log as a local diffeomorphism. bounds ensure JGblog (x) remains invertible on Rϵ (b p), establishing G b log and prove that We are then able to apply analogous arguments for the proof Theorem 6 to G log log m b b Φ := limm→∞ (G ) exists and is surjective from ∂Rϵ (b p) onto R(b p). ■ Proof of Theorem 8. We use the same definition of N3 and En (c) in the proof of Lemma 7, where the event En (c) occurs with probability at least 1 − n−c when n ≥ N3 . Assume En (c) holds. Then similar to Theorem 4, we have dH (R(b p), R(p)) ≤

2ϵlog ξ (n) γlog

.

(110) 2

As shown in the proof of Lemma 7, the adaptive step size αn satisfies inf x∈D αn (x) ≥ CK2h . Similar to Theorem 5 and Corollary 3, we can show that 8ϵ m log dH (Xbm , R(b p)) ≤ ρ , γlog n

(111)

C γ

where ρn = 1 − K16log h2 By the triangle inequality for the Hausdorff metric, log log dH (Xbm , R(p)) ≤ dH (Xbm , R(b p)) + dH (R(b p), R(p)).

Substituting (110) and (111) into this inequality yields the total error bound.

40

(112) ■

Related documents

• The Iterates of Nesterov's Accelerated Algorithm Converge in The Critical Regimes#1108962
• A trust region method for minimization of nonsmooth functions with linear constraints#429984
• Rigidity of reflectionless Riesz measures and Sobolev estimates for Riesz transforms on Lipschitz graphs#1116510
• Lyapunov Constructions for System Interconnections Arising from Adaptation in Some Optimization Methods#625400
• Lyapunov Constructions for System Interconnections Arising from Adaptation in Some Optimization Methods#623127