Riemannian Simultaneous Inference for Tangent Vector Field Regression
arXiv:2609.21910v1 [stat.ML] 18 Sep 2026
Xiaotian Chang1 1 2
Yangdi Jiang1
Qirui Hu2∗
Division of Mathematical Sciences, Nanyang Technological University
School of Statistics and Data Science, Shanghai University of Finance and Economics
Abstract We consider nonparametric tangent vector field regression on a Riemannian manifold without boundary. Because responses at different points lie in different tangent spaces, the proposed kernel estimator first parallel transports nearby responses to the target tangent space and then forms a volume-corrected local average. We first derive its uniform second-order bias, finitebandwidth covariance, and stochastic rate. For simultaneous inference, the tangent norm is written as a supremum over the unit tangent bundle. Exact covariance whitening gives a unitvariance Gaussian field whose correlation length is of order h along the base manifold and of order one along the fibre. Its local covariance geometry leads to a Gumbel limit with an explicit intrinsic constant. Combining this limit with Gaussian approximation and cross-fitted covariance estimation yields a feasible simultaneous confidence tube for the regression field. We further discuss improved finite-sample inference with bandwidth selection and high-order bias corrections. Simulations on various manifolds support the proposed inference procedure. A randomized reconstruction of global wind data illustrates how the tube’s cross-sections describe spatially varying uncertainty.
Keywords: Gaussian extremes; kernel regression; parallel transport; Riemannian manifolds; simultaneous confidence tubes. ∗
Corresponding author: [email protected].
1
1
Introduction
Many scientific measurements describe local motion on a curved domain. Their locations lie on a Riemannian manifold M, while their values lie in the tangent spaces at those locations. Horizontal wind provides a familiar example: a geographic location is a point x ∈ S 2 , and the wind velocity is an element of the tangent plane Tx S 2 . Estimating the mean circulation and quantifying uncertainty over the globe therefore requires statistical methods for an entire tangent vector field, rather than separate analyses of vectors in one fixed Euclidean space (Fan et al., 2018; Robert-Nicoud et al., 2024). The same data structure occurs in robotics. The full orientation of an end effector is a rotation R ∈ SO (3), and a body angular velocity ω ∈ R3 determines the tangent velocity R [ω]× ∈ TR SO (3), where [ω]× v = ω × v. Orientation-dependent dynamical systems are used to describe end-effector motion (Beik Mohammadi et al., 2024). Repeated orientation–velocity observations can therefore be used to estimate mean rotational motion and its uncertainty as orientation varies. Figure 1 illustrates these two motivating settings1 ; the robot panel is schematic. We further present the data analysis concerning planetary-scale wind reconstruction in Section 6. To place these applications in a common statistical framework, we observe independent copies of a pair (X, Y ) satisfying
Xi ∈ M,
Yi ∈ TXi M,
Yi = V (Xi ) + ηi ,
E (ηi | Xi ) = 0,
(1.1)
where X has density p with respect to Riemannian volume, and the unknown regression field is V (x) = E (Y | X = x) ∈ Tx M. Thus V is a section of the tangent bundle: it assigns a vector to each location in its own tangent space. The conditional error covariance may vary with location. Our first objective is to estimate V uniformly over M using a kernel estimator. Our principal objective is simultaneous inference: given α ∈ (0, 1), construct regions Cn,1−α (x) ⊂ Tx M such that
P {V (x) ∈ Cn,1−α (x) for every x ∈ M} −→ 1 − α. 1 The wind image uses NCEP–NCAR Reanalysis 1 on a 2.5◦ grid, provided by NOAA PSL, Boulder, Colorado, USA (https://psl.noaa.gov); land imagery is from Natural Earth. This snapshot differs from the data in Section 6.
2
Orientation R ∈ SO (3); body angular velocity ω ∈ R3 . Tangent response: Y = R [ω]× ∈ TR SO (3). Target: mean tangent velocity V (R) = E (Y | R), with simultaneous uncertainty.
Figure 1: Tangent-vector observations on a sphere and a rotation group. Left: 500-hPa wind on 1 July 2024, 12:00 UTC; color indicates speed and arrows indicate direction. Right: schematic end-effector orientations and angular velocities on SO (3). These regions form a confidence tube around the estimated field, covering all locations and tangent directions simultaneously. We allow a nonparametric mean field and spatially varying error covariance on a known compact manifold. Kernel methods on manifolds provide a natural starting point. Volume-corrected kernel density estimation and scalar-response regression on closed manifolds were studied by Pelletier (2005, 2006); Henry and Rodriguez (2009) developed robust regression in the same geometric setting. Tangent responses introduce a further issue: both the location and the space containing the response vary. Singer and Wu (2012) use local alignments to construct vector diffusion maps and relate their limit to the connection Laplacian, with applications to vector-field interpolation and regression. Their manifold-learning setting differs from the known-manifold random-design model considered here. Das and Snasel (2026) study concentration for bundle-valued statistics transported to a reference fibre. Such concentration results address a different objective from calibrating a shrinkingbandwidth regression maximum over all target fibres. A complementary literature models tangent fields through Gaussian processes. On the sphere, Fan et al. (2018) construct covariance models from surface gradients and curls of scalar potentials,
3
with likelihood-based parameter estimation. Hutchinson et al. (2021) develop gauge-independent projected kernels, whereas Robert-Nicoud et al. (2024) construct intrinsic Hodge–Matérn Gaussian vector fields. The latter two constructions provide Gaussian process priors for Bayesian prediction and uncertainty quantification: the unknown field is modeled as random, and inference uses its conditional distribution given the observations. Our objective is instead frequentist simultaneous coverage for a fixed unknown regression field under independent noisy observations; the Gaussian field below is a reference law for estimation error, not a prior on V . Simultaneous inference has a long history in Euclidean statistics. Bickel and Rosenblatt (1973) study global deviations of density estimators, and Johnston (1982) obtain maximal-deviation limits for nonparametric regression estimators. Härdle and Marron (1991) construct simultaneous regression error bars by residual bootstrap, providing an alternative to analytic extreme-value calibration. Modern Gaussian approximation and anti-concentration methods support confidence-band construction without requiring an extreme-value limit (Chernozhukov et al., 2014, 2015, 2022). For manifold-indexed Gaussian fields, Qiao and Polonik (2018) study extrema under rescaling, and Qiao (2021a) develops excursion and extreme-value results for locally stationary Gaussian and chi fields. Related statistical uses include confidence regions for density ridges (Qiao, 2021b) and multivariate random-field inference (Taylor and Worsley, 2008). Two difficulties remain in the present problem. First, nearby responses must be brought into a common tangent space before averaging. A global orthonormal frame need not exist, and coordinatewise confidence bands would depend on the frames chosen. Curvature and a nonuniform design also affect the local mean and covariance. Second, simultaneous inference requires the dependence between estimation errors at different locations. Gaussian approximations defined separately in local frames must agree where those frames overlap. Their excursion probabilities must then be combined to obtain one threshold for all locations and tangent directions. We address the first difficulty by transporting nearby responses along short geodesics before taking a volume-corrected kernel average. Its second-order bias identifies the effects of the covariant derivatives of V and the design density. Uniform stochastic bounds then separate the local noise, transported-signal variation and denominator error. For simultaneous inference, we construct one Gaussian reference section whose covariance matches the leading estimation error at every pair of locations. A shared Hilbert-space repre4
sentation ensures that its local representations agree on chart overlaps. Expressing its norm as a maximum over unit tangent directions allows us to compare its maximum with the regression maximum on a common finite grid. A local covariance expansion then determines high-excursion probabilities, which combine into a global Gumbel law with an explicit constant. This gives an analytic critical value for the confidence tube. To make the confidence tube feasible, we estimate the error covariance from cross-fitted residuals at each observation location. These covariance estimates are transported to the target tangent space and combined with squared kernel weights, as required for the covariance of a weighted mean. Separate bandwidths for regression, residual prediction and covariance smoothing make the covariance estimate sufficiently accurate for this calibration. Undersmoothing gives coverage for V ; a separate target-local construction gives inference for the smoothed field when spatial resolution is part of the estimand. Simulations and the wind reconstruction illustrate the resulting confidence tubes. The rest of the paper is organized as follows. Section 2 introduces the model, geometry, estimator and assumptions. Section 3 establishes the estimation results, and Section 4 develops simultaneous inference. Sections 5 and 6 present the simulations and wind reconstruction. The supplementary material follows the main results in proof order and provides implementation details.
2
Methodology
We first specify the geometry needed to compare responses at different locations and then define the kernel estimator. The geometry is assumed known. Let (M, g) be a compact, connected, smooth Riemannian manifold of dimension d ≥ 2 without boundary. Write ρ for geodesic distance, Vol for Riemannian volume, ⟨·, ·⟩x and ∥·∥x for the inner product and norm on Tx M, and ιM > 0 for the injectivity radius. The exponential map expx sends an initial tangent vector u to the time-one endpoint of the geodesic starting at x with velocity u. On a tangent ball of radius less than ιM , it is a diffeomorphism onto a normal neighborhood of x. Its inverse is the logarithm map logx . If ρ (x, q) < ιM , let Pq→x denote Levi–Civita parallel −1 = P transport along the unique minimizing geodesic. Transport is an isometry and Pq→x x→q .
5
For q = expx (u) in a normal neighborhood, define the volume density by
d Vol (q) = Θx (q) du,
where du is Lebesgue measure induced by the tangent inner product. On uniformly small normal neighborhoods, Θx (q) is bounded above and away from zero. Its reciprocal removes volume distortion in the kernel calculation. We write ∇ for the Levi–Civita covariant derivative and Γ (T M) for the sections of the tangent bundle. The endomorphism bundle End (T M) has at x the fibre of linear maps from Tx M to itself; a section assigns one such map to every point. The covariance field below is an example. A local orthonormal frame supplies coordinates, but the n
o
section and its norm do not depend on that frame. For example, Tx S 2 = u ∈ R3 : x⊤ u = 0 and n
o
TR SO (3) = RA : A⊤ = −A . The unit tangent bundle, needed later for simultaneous inference, is
S (T M) = {(x, v) : x ∈ M, v ∈ Tx M, ∥v∥x = 1} . Its base has dimension d and its fibres have dimension d − 1. Only local product coordinates are used; a global frame is not assumed. See Lee (2018) for these geometric conventions. Under model (1.1), neither the density p nor the regression field V is known. The conditional error covariance is the self-adjoint operator
Σ (x) = E {ηi ⊗ ηi | Xi = x} : Tx M → Tx M,
(u ⊗ v) w = ⟨v, w⟩x u.
To estimate at x, we first restrict attention to observations within geodesic distance h from x, transport their responses into Tx M, and then average them there. The factor Θx (q)−1 removes the normal-coordinate volume distortion from the local mean. This is the geometric counterpart of a Nadaraya–Watson average. Let K : [0, ∞) → R be radial and supported on [0, 1], and write ⋆
K (u) = K (∥u∥) ,
Z
R (K) =
⋆
2
{K (u)} du,
Rd
6
1 µ2 (K) = d
Z Rd
∥u∥2 K ⋆ (u) du.
(2.1)
For 0 < h < ιM , define
Wh (x, q) =
Θx (q)−1 K {ρ (x, q) /h} ,
ρ (x, q) < ιM ,
0,
ρ (x, q) ≥ ιM .
Because K is supported on [0, 1], Wh (x, q) = 0 unless ρ (x, q) ≤ h, where the short-geodesic transport is uniquely defined. The empirical denominator and transported numerator are
pbh (x) =
n 1 X Wh (x, Xi ) , nhd i=1
bh (x) = N
n 1 X Wh (x, Xi ) PXi →x Yi ∈ Tx M. nhd i=1
Let pn > 0 be a deterministic sequence with pn ↓ 0. For measurability on every sample, put pb†h (x) = pbh (x) ∨ pn and define Vbh (x) =
bh (x) N
pb†h (x)
∈ Tx M.
When the denominator floor is inactive, the normalized weights Wh (x, Xi ) /
(2.2) P
j Wh (x, Xj ) sum to
one. Under the assumptions below, the floor is uniformly inactive with probability tending to one. To analyze this local average, we introduce its population counterpart:
ph (x) = E {pbh (x)} , 1 Wh (x, q) Pq→x V (q) p (q) d Vol (q) , hd M Nh (x) Vh (x) = . ph (x) Z
Nh (x) =
The smoothed field Vh is the population target of the local average and separates stochastic error from smoothing bias. The first three assumptions below justify the local mean, covariance and Gaussian calculations; the fourth makes Vh − V negligible at the Gumbel scale for simultaneous inference on V . Conditions for the separate pilot and covariance bandwidths appear in Section 4.3. We first present and discuss some necessary assumptions for subsequent theoretical analysis.
Assumption 2.1. The function K ⋆ (u) = K (∥u∥) is nonnegative, belongs to Cc4 Rd , is supported 7
on the closed unit ball, and satisfies
R
Rd K
⋆ (u) du = 1, 0 < R (K) < ∞.
Assumption 2.2. The design distribution has density p ∈ C 3 (M) with respect to Vol, and there are constants 0 < pmin ≤ pmax < ∞ such that pmin ≤ p (x) ≤ pmax , x ∈ M. The regression field satisfies V ∈ C 3 (Γ (T M)). Assumption 2.3. Choose an everywhere-defined Borel version of the conditional law of η given X = x. Its conditional mean is zero, and its covariance field Σ is a self-adjoint C 3 section of End (T M). There are constants Cη < ∞ and 0 < λmin ≤ λmax < ∞ such that, for every x ∈ M, every unit a ∈ Tx M, and every t ∈ R,
E {exp (t ⟨a, η⟩x ) | X = x} ≤ exp Cη2 t2 /2 ,
and λmin ∥u∥2x ≤ ⟨u, Σ (x) u⟩x ≤ λmax ∥u∥2x ,
u ∈ Tx M.
Assumption 2.4. The deterministic bandwidth h = hn satisfies
h ↓ 0,
(log n)5 −→ 0, nhd
nhd+4 log (1/h) −→ 0.
(2.3)
The smoothness and tail conditions support uniform inference. In the scalar-response setting, Pelletier (2006) uses a closed manifold, a normalized compactly supported kernel, positive design density, and twice continuously differentiable density and regression functions, with bounded responses for the pointwise expansions. Our C 3 fields and C 4 kernel provide the additional derivative control needed for covariance whitening and moving kernel overlaps. The sub-Gaussian assumption permits unbounded Gaussian errors but imposes uniform tail control for the high-dimensional maximum comparison (Chernozhukov et al., 2022); ellipticity makes that whitening stable. Local covariance nondegeneracy and regularity are also central in Gaussian extreme-value theory (Qiao and Polonik, 2018; Qiao, 2021a). These conditions allow smooth spatially varying anisotropic co
variance. For example, the normalized nonnegative kernel K ⋆ (u) = cd 1 − ∥u∥2
5
1 {∥u∥ ≤ 1}
satisfies the kernel condition and can be combined with smooth tangent fields and conditionally Gaussian errors having smooth uniformly positive covariance.
8
The last condition in (2.3) undersmooths the second-order bias so that the confidence tube is centered at V . For h = n−α , the two rate conditions hold when 1 1 <α< . d+4 d
(2.4)
Similar undersmoothing conditions occur in scalar simultaneous regression inference. For example, Assumption (A5) of Cai et al. (2021) allows h = n−α with 1/5 < α < 1/3, a subrange of (2.4) when d = 1. For scalar regression on Riemannian manifolds, Henry and Rodriguez (2009) use nhd → ∞ and n1/(d+4) h → 0 for a centered pointwise normal limit (Assumption A3 with β = 0 and Theorem 4.1). For h = n−α , these bandwidth conditions give the same range (2.4). On a flat Euclidean domain, parallel transport is the identity, Θx (q) = 1, and (2.2) is the usual vector-valued Nadaraya–Watson estimator. The formulas below then reduce to the familiar coordinatewise bias and covariance expressions.
3
Estimation theory
The estimator admits an exact decomposition into smoothing bias, random-design fluctuation, and transported noise. We use it to derive the uniform bias and covariance expansions required for simultaneous inference. Define the centered transported-signal and noise processes in Tx M by n 1 X Rn,h (x) = √ Wh (x, Xi ) {PXi →x V (Xi ) − Vh (x)} , nhd i=1 n 1 X Un,h (x) = √ Wh (x, Xi ) PXi →x ηi . nhd i=1
(3.1)
Then, on every sample, √
√ n
o
nhd Vbh (x) − Vh (x) =
Rn,h (x) + Un,h (x) +
n
o
nhd pbh (x) − pb†h (x) Vh (x)
pb†h (x)
.
(3.2)
Here Rn,h is the fluctuation of the transported signal under random design, and Un,h is the leading noise term. The final numerator term accounts for denominator truncation. For all sufficiently
9
large n, it vanishes on Ep,n = {inf x∈M pbh (x) ≥ pmin /2}, whose probability tends to one. We first calculate the smoothing bias Vh − V in normal coordinates, then the covariance of Un,h . Empirical-process bounds over a finite atlas control the random terms uniformly over M. For a scalar function f , ∇f is its Riemannian gradient and ∆f = div(∇f ) is its Laplace– Beltrami operator. For vector fields, ∇U V denotes the covariant derivative of V along U ; in particular, ∇∇p V differentiates V along the density gradient. The connection Laplacian is the trace of the second covariant derivative: ∆∇ V (x) =
d n X
o
∇ej ∇ej V − ∇∇ej ej V (x) ,
j=1
where {e1 , . . . , ed } is any local orthonormal frame. Lemma 3.1. Under Assumptions 2.1–2.2, uniformly in x ∈ M, h2 µ2 (K) ∆p (x) + O h3 , 2 1 1 ∇ 2 ∆ V (x) + ∇∇p(x) V (x) + O h3 . Vh (x) − V (x) = h µ2 (K) 2 p (x)
ph (x) = p (x) +
Under Assumption 2.4,
sup |pbh (x) − p (x) | = OP
x∈M
h2 +
q
s
log (1/h) , nhd
sup ∥Rn,h (x)∥x = OP h log (1/h) .
x∈M
The leading random fluctuation is determined by the covariance of Un,h in (3.1):
Ch (x) =
i 1 h 2 E W (x, X) P Σ (X) P . X→x x→X h hd
Its conditional-design analogue, still involving the unknown Σ, is
Cbh (x) =
n 1 X Wh (x, Xi )2 PXi →x Σ (Xi ) Px→Xi . nhd i=1
By construction, Cbh (x) = Cov {Un,h (x) | X1 , . . . , Xn } .
10
(3.3)
The finite-bandwidth population covariance scale for the leading regression noise and its limit are Ωh (x) =
Ch (x) , ph (x)2
Ω (x) =
R (K) Σ (x) . p (x)
Lemma 3.2. Under Assumptions 2.1–2.3, uniformly in x,
sup Cbh (x) − Ch (x) x∈M
op
= OP
s log (1/h)
nhd
log (1/h) + . nhd
Moreover,
Ch (x) = p (x) R (K) Σ (x) + O h2 ,
(3.4)
and sup Ω (x)−1/2 {Ωh (x) − Ω (x)} Ω (x)−1/2 x∈M
op
= O h2 .
(3.5)
The volume correction removes the volume-density term from the leading bias, leaving the connection Laplacian and the interaction between the design density and the field. The covariance expansion instead retains one inverse volume factor because the kernel weight is squared. Its leading term R (K) Σ/p determines the local shape of the confidence tube. Combining these expansions with (3.2) gives the uniform rate. Theorem 3.1 (Uniform rate). Under Assumptions 2.1–2.4,
sup Vbh (x) − V (x) x∈M
x
= OP
s
h2 +
log (1/h)
nhd
.
The stochastic term has the usual effective-sample-size scale nhd , with a logarithmic factor for uniformity over the manifold. For inference, this rate alone is not sufficient: the maximum fluctuates on the smaller Gumbel scale. Section 4 therefore retains the leading noise field and controls the remaining terms separately at that scale. In a flat chart, away from a boundary, the leading bias is
h2 µ2 (K)
1
2
∆V +
d X j=1
11
(∂j log p) ∂j V
,
and the covariance is R (K) Σ/p: these are the vector-valued Nadaraya–Watson expressions. The same formulas hold globally on a flat torus. Thus curvature does not change the powers h2 and
nhd
−1/2
; it enters transport, the volume correction and higher-order terms. This agrees with
the Euclidean-order bias, variance and integrated rates for scalar manifold regression in Pelletier (2006).
4
Simultaneous inference
Simultaneous inference requires a critical value for the largest standardized estimation error. We obtain it in three steps: approximate the regression maximum by a Gaussian maximum, derive its Gumbel law, and show that replacing the unknown covariance by an estimate preserves that law.
4.1
Gaussian reference field and distributional comparison
To analyze the maximum, we express a vector norm as a scalar supremum: for every tangent field G, sup ∥G (x)∥x =
x∈M
sup x∈M, v∈Tx M ∥v∥x =1
⟨v, G (x)⟩x .
(4.1)
The unit tangent bundle S (T M) was defined in Section 2. Equation (4.1) converts the vector-norm problem to a scalar supremum without choosing a global frame. Apply (4.1) to the standardized noise field G (x) = Ch (x)−1/2 Un,h (x), where Un,h and Ch are defined in (3.1) and (3.3). Its directional coordinate is D
Zn,h (x, v) = v, Ch (x)−1/2 Un,h (x)
E x
,
(x, v) ∈ S (T M) .
By (3.1), this coordinate is a weighted sum of transported errors. Its covariance can be written as an L2 inner product of the corresponding kernel, transport, and covariance features. We therefore place all such features in one Hilbert space and apply a single isonormal Gaussian process. This construction produces a scalar Gaussian field on S (T M) whose finite-bandwidth covariance agrees exactly with that of the whitened empirical noise. All quantities are defined intrinsically on the tangent bundle. Let H
= L2 {M, T M, p d Vol} be the Hilbert space of square-integrable measurable 12
tangent sections. R
An element is a field f assigning f (q) ∈ Tq M to almost every q, with
2 M ∥f (q)∥q p (q) d Vol (q) < ∞; fields equal almost everywhere under p d Vol represent the same
element. The inner product integrates the tangent-space inner products:
⟨f, g⟩H =
Z M
⟨f (q) , g (q)⟩q p (q) d Vol (q) .
For x ∈ M, define the feature operator Fh,x : Tx M → H by (Fh,x v) (q) = h−d/2 Wh (x, q) Σ (q)1/2 Px→q Ch (x)−1/2 v.
(4.2)
The factors in (4.2) have distinct roles. Starting with a direction v at x, the operator Ch (x)−1/2 standardizes the noise variance, parallel transport moves that direction to a possible observation location q, and Σ (q)1/2 accounts for the error covariance there. The kernel weight specifies how strongly that observation contributes. Inner products of these features reproduce the covariance between standardized directional errors at different locations. ∗ : H → T M denotes the Hilbert-space adjoint, defined by Here Fh,x x
D
∗ ⟨Fh,x v, f ⟩H = v, Fh,x f
E x
,
v ∈ Tx M,
f ∈H.
∗ F Using the definition of Ch gives ⟨Fh,x v, Fh,x w⟩H = ⟨v, w⟩x . Equivalently, Fh,x h,x = ITx M , so the
feature map is an exact isometry. Let W be an isonormal Gaussian process over H : this is a centered jointly Gaussian family indexed by f ∈ H , linear in f , with E {W (f ) W (g)} = ⟨f, g⟩H . A single such process supplies the randomness at every location. Define
Zh (x, v) = W (Fh,x v) ,
(x, v) ∈ S (T M) .
This field has variance one and exactly the covariance of the empirical coordinate Zn,h defined above. The scalar field also defines a Gaussian random section of T M. For a local orthonormal frame
13
e1 (x) , . . . , ed (x), set
Zh (x) =
d X
W {Fh,x ej (x)} ej (x) ,
Zh (x, v) = ⟨v, Zh (x)⟩x .
(4.3)
j=1
An orthogonal change of local frame changes the Gaussian coordinates and basis vectors by inverse transformations, leaving the sum unchanged. The local definitions therefore agree on chart overlaps. Feature regularity gives a continuous version, so Zh is a random element of the space of continuous sections. It is indexed by x ∈ M and takes values in Tx M; Zh is its scalar representation indexed by (x, v) ∈ S (T M). We next quantify how the estimator inherits the Gaussian calibration. The relevant normalization is ah =
q
2d log (1/h)
(4.4)
because the Gaussian maximum is of order ah , while its fluctuations are of order a−1 h . Accordingly,
an additive perturbation must be oP a−1 , and a relative covariance perturbation must be oP a−2 . h h The Gaussian comparison itself is distributional. For the empirical coordinate Zn,h , independence and conditional centering of the observations give the exact covariance identity
Cov {Zn,h (x, v) , Zn,h (y, w)} = ⟨Fh,x v, Fh,y w⟩H = Cov {Zh (x, v) , Zh (y, w)} . ∗ F In particular, the cross-covariance operator of the Gaussian section is Fh,x h,y : Ty M → Tx M; at
x = y it is the identity. Consequently Ωh (x)1/2 Zh (x) has the covariance of ph (x)−1 Un,h (x), the leading scaled regression noise. This covariance match follows from the sampling model and does not require Gaussian observations. The remaining issue is whether it also yields an approximation to the distribution of the maximum. The target-centered oracle regression maximum is V Tn,h = sup Ω (x)−1/2
√
n
nhd Vbh (x) − V (x)
x∈M
o x
.
We compare both fields on the same deterministic grid and control the error from replacing the
14
−3 continuum by that grid. In a finite bundle atlas, use base mesh ha−3 h and fibre mesh ah . The 3(2d−1)
grid Gh has |Gh | ≲ h−d ah
points. Uniform scaled first derivatives of both fields are OP (ah ),
, smaller than the fluctuation scale a−1 so interpolation changes either maximum by OP a−2 h . h At each grid point, the empirical field is a normalized sum of independent centered scores. Their
variance is one, their sub-exponential norm is O h−d/2 , and their fourth moment is O h−d . Theorem 2.1 of Chernozhukov et al. (2022) therefore gives
sup P max Zn,h ≤ t − P max Zh ≤ t t
Gh
Gh
"
log5 {n|Gh |} ≤C nhd
#1/4
= o (1) .
Gaussian anti-concentration (Chernozhukov et al., 2015) controls the threshold shifts needed to −3/2
return to the continuum. For example, a shift of ah
dominates either interpolation error with −1/2
probability tending to one and changes the Gaussian grid cdf by O ah
.
Finally, the exact decomposition in Section 3 returns us from the noise maximum to the regression estimator. With Lh = log (1/h), V ah Tn,h − sup Zn,h = OP Lh S(T M)
s
h2 +
Lh Lh + + hLh + nhd nhd
q
nhd+4 Lh = oP (1) .
The final term is the smoothing bias. The other terms arise from normalization and transportedsignal variation under random design. Anti-concentration also transfers this small estimator perturbation to a cdf comparison. The supplementary material supplies these arguments, yielding the following proposition. Proposition 4.1. Under Assumptions 2.1–2.4, (
sup P t∈R
)
(
sup Zn,h ≤ t − P
S(T M)
)
sup Zh ≤ t
−→ 0,
S(T M)
and V ah Tn,h − sup Zn,h = oP (1) . S(T M)
Consequently, sup P t∈R
n
o
V Tn,h ≤t
(
−P
)
sup Zh ≤ t S(T M)
15
−→ 0.
It remains to find a critical value for the Gaussian maximum. The next subsection derives its Gumbel limit from the local covariance structure; Proposition 4.1 then transfers this calibration to the regression maximum.
4.2
Local covariance and Gaussian calibration
The Gaussian maximum depends on how quickly nearby indices decorrelate, not just on their unit marginal variances. Moving a location changes the kernel neighborhood; rotating a tangent direction changes which component of the vector is measured. Both changes enter the same Taylor expansion. For rh {(x, v) , (y, w)} = Cov {Zh (x, v) , Zh (y, w)}, exact whitening gives
1 − rh {(x, v) , (y, w)} =
1 ∥Fh,x v − Fh,y w∥2H . 2
(4.5)
Thus the quadratic covariance expansion is obtained by differentiating the feature map and taking its Hilbert-space Gram matrix. The two kinds of displacement have different scales. A spatial displacement of size h changes the translated kernel by order one, whereas the fibre sphere is not shrinking. Using R (K) from (2.1), define R
cK =
⋆ 2 Rd (∂1 K ) du
2R (K)
> 0.
(4.6)
At a center s = (x, v), choose orthonormal base and fibre coordinates and write the bundle chart directly as Φh,s (hu, z), where u ∈ Rd and z ∈ Rd−1 . Here hu is the physical base coordinate. The fibre coordinate is adjusted linearly as the base point moves so that the two feature-derivative blocks are orthogonal at s. Under Assumptions 2.1–2.3, the joint expansion is 1 − rh Φh,s (hu, z) , Φh,s hu′ , z ′
=
u − u′
z − z′
⊤
cK Id + Eh,s
0
u − u′
0 1 2 Id−1
16
z − z′
′ ′ + Rh,s u, z, u , z ,
(4.7)
where sups ∥Eh,s ∥op ≤ Ch and, on a fixed sufficiently small rescaled coordinate neighborhood, |Rh,s u, z, u′ , z ′ | ≤ C ∥u∥ + u′ + ∥z∥ + z ′
n
u − u′
2
+ z − z′
2
o
.
The constants and neighborhood are uniform in s and sufficiently small h. The spatial derivative block includes a factor h by the chain rule; the fibre derivative block does not. The factor 1/2 multiplying the Hilbert squared distance accounts for the fibre block Id−1 /2. The mixed block is exactly zero at the chart center after the adjustment. It need not vanish away from the center: its variation is included in the uniform remainder. For excursion probabilities, we need this expansion uniformly around every nearby reference point, not only at the chart center. For s = (x, v) and rescaled coordinates t = (tu , tz ), define the (2d − 1) × (2d − 1) matrix directly from the correlation function:
[Qh,s (t)]ab =
1 ∂2 rh Φh,s (htu , tz ) , Φh,s ht′u , t′z ′ 2 ∂ta ∂tb
. t′ =t
Equivalently, this is one half of the Gram matrix of the derivatives of Fh,x v in these coordinates. For a small coordinate displacement ∆, the quadratic form ∆⊤ Qh,s (t) ∆ gives the leading loss of correlation. A larger value means that the Gaussian field decorrelates more quickly in that direction. The following statement gives the uniform control needed to use this calculation throughout the bundle. Lemma 4.1. Under Assumptions 2.1–2.3, there are constants r0 , c, C > 0, independent of s and sufficiently small h, such that the adjusted charts above are defined on a common rescaled ball Br0 ⊂ R2d−1 . For t, t′ ∈ Br0 ,
cI2d−1 ⪯ Qh,s (t) ⪯ CI2d−1 ,
Qh,s (t) − Qh,s t′
op ≤ C
t − t′ .
(4.8)
At t = 0, Qh,s (0) is the block matrix in (4.7). For any a, b, t ∈ Br0 , 1 − rh {Φh,s (hau , az ) , Φh,s (hbu , bz )} − (a − b)⊤ Qh,s (t) (a − b) ≤ C max {∥a − t∥ , ∥b − t∥} ∥a − b∥2 .
17
(4.9)
The lower bound in (4.8) ensures that every base and fibre direction contributes to the covariance loss. The Lipschitz bound and (4.9) allow the exact matrix at any reference point to replace nearby matrices with a uniformly smaller-order error. At the center, replacing cK Id + Eh,s by cK Id adds at most Ch ∥au − bu ∥2 . The supplementary material constructs the adjustment and proves Lemma 4.1. Compact kernel support gives rh {(x, v) , (y, w)} = 0 when ρ (x, y) > 2h. Away from the rescaled diagonal there is also a uniform one-sided gap below one. The qualification “one-sided” matters because opposite directions in the same fibre have correlation −1. The exact covariance depends on p and Σ, but whitening removes them from the leading quadratic blocks. This explains why the leading excursion constant depends on the kernel, dimension and Riemannian volume. To turn the covariance expansion into an extreme-value law, we use the local-stationarity ideas developed for manifold-indexed fields by Qiao and Polonik (2018); Qiao (2021a), with the spatial and fibre scales verified here on the unit tangent bundle. We first calculate the probability of a high excursion within one small spatial neighborhood, allowing all tangent directions there. The local coefficient is
q
det Qh,s (t): rapid decorrelation in
more directions increases the contribution to the excursion probability. To combine these coefficients over a region A ⊂ S (T M), partition A, up to chart boundaries, into disjoint pieces represented by coordinate domains Dν in the rescaled charts centered at sν , and set d
Ih (A) = h
XZ ν
q
Dν
det Qh,sν (t) dt.
The determinant transforms with the coordinate Jacobian, so the sum does not depend on the chosen partition or charts. The factor hd compensates for rescaling the base coordinates. Thus Ih records the integrated local excursion coefficient using the covariance expansion already obtained. Partition M into mh ≍ h−d regular cells Jk,h of diameter comparable to h, and let π (x, v) = x be the bundle projection. For a fixed small δ > 0, shrink each cell by a fraction δ in its simplex
δ = π −1 J δ coordinates and lift it to Ek,h k,h , including the full fibre sphere. The removed boundary
strips are the collars. Let Dhδ consist of these lifted shrunken cells and the lifted pieces of a regular triangulation of their collars. The supplementary material gives a construction with uniform shape and boundary bounds after rescaling base coordinates by h−1 . Proposition 4.2. Under Assumptions 2.1–2.3, fix a sufficiently small δ > 0. As u → ∞, uniformly 18
over Ah ∈ Dhδ and small h, (
)
Zh (x, v) > u
sup
P
= {1 + o (1)} π −(2d−1)/2 u2d−1 Ψ (u) h−d Ih (Ah ) ,
(4.10)
(x,v)∈Ah
where Ψ (u) = u−1 (2π)−1/2 e−u /2 . 2
Proposition 4.2 turns the local covariance geometry into a cell exceedance probability. Its uniform relative error is essential because the number of cells diverges as h ↓ 0. The proof in the supplementary material controls chart boundaries and double excursions, using the uniform bounds of Dębicki et al. (2017). The center blocks in (4.7) also give
d/2
Ih {S (T M)} −→ cK 2−(d−1)/2 Vol (M) Vol S d−1 .
Consequently, summing (4.10) over the manifold produces an intensity proportional to 2
h−d u2d−2 e−u /2 . With ah as in (4.4), put Cvec =
2d/2 dd−1 d/2 c Vol (M) . π d/2 Γ (d/2) K
(4.11)
The constant Cvec combines the base-volume and fibre-area integrals above. The classical expanded Gumbel centering is βhvec = ah +
(d − 1) log log (1/h) + log Cvec . ah
(4.12)
To see why this intensity gives a Gumbel law, fix z ∈ R and set uh,z = βhvec + z/ah . Count the shrunken cells that contain an exceedance:
Wh,δ =
mh X k=1
1
sup Zh > uh,z
E δ
,
λh,δ = EWh,δ .
k,h
Summing the uniform cell probabilities gives the mean of this count. Cells more than 2h apart have independent Gaussian fields, so each cell depends on only a bounded number of neighbors. The correlation gap makes joint exceedances in neighboring shrunken cells negligible. The Chen–Stein
19
bound of Arratia et al. (1989) then yields, for fixed δ and z as h ↓ 0, λh,δ = e−z + O (δ) + o (1) ,
(4.13)
dTV {L (Wh,δ ) , Po (λh,δ )} −→ 0. Here dTV is total variation distance and Po (λ) denotes the Poisson law with mean λ. The probability of an exceedance in the omitted collars is O (δ) + o (1). Thus (4.13) converts the event of no exceedance into (
P
)
sup Zh ≤ uh,z
= exp −e−z + O (δ) + o (1) .
S(T M)
Letting h ↓ 0 first and then δ ↓ 0 proves the following theorem. The supplementary material supplies the joint-tail bounds, Poisson approximation and collar calculation. Theorem 4.1 (Gaussian sphere-bundle Gumbel law). Under Assumptions 2.1–2.3, as h ↓ 0 with h < ιM /2, for every z ∈ R, "
P ah
(
)
sup (x,v)∈S(T M)
Zh (x, v) − βhvec
#
≤ z −→ exp −e−z .
(4.14)
The full fibre sphere contains both scalar signs, and its area therefore supplies the intensity for the tangent norm. Combining Proposition 4.1 with Theorem 4.1 gives, for every z ∈ R, h
n
o
i
V P ah Tn,h − βhvec ≤ z −→ exp{−e−z }.
On a flat torus, this is the joint norm limit for periodic Rd -valued kernel regression. The corresponding Euclidean result holds on a compact inference region inside the design support. Proposition 4.3 (Euclidean specialization). Let d ≥ 1, D =
Qd
j=1 [ℓj , rj ] with ℓj < rj , and consider
independent copies of (X, Y ) satisfying Y = V (X) + η, where X, Y ∈ Rd . Suppose Assumptions 2.1 and 2.4 hold, and the Euclidean versions of Assumptions 2.2–2.3 hold uniformly on a fixed open neighborhood of D. Use the Euclidean estimator (2.2), retaining observations outside D. With |D|
20
denoting Lebesgue volume, set √ V Tn,h,D = sup Ω(x)−1/2 nhd {Vbh (x) − V (x)} ,
CD =
x∈D
βh,D = ah +
(d − 1) log log(1/h) + log CD , ah
ah =
q
2d/2 dd−1 d/2 c |D|, π d/2 Γ(d/2) K
2d log(1/h).
Then, for every z ∈ R, h
i
V P ah {Tn,h,D − βh,D } ≤ z −→ exp{−e−z }.
(4.15)
The limit also holds with Ω replaced by a continuous positive-definite estimate whose uniform relative error is oP (a−2 h ), as in (4.22) with the supremum taken over D. For d = 1, write D = [ℓ, r], L = r − ℓ and Σ(x) = σ 2 (x). Then p
nh p(x) b V Tn,h,D = sup p |Vh (x) − V (x)|, R(K) σ(x) x∈D
√ L 2cK CD = . π
+ There is no log log(1/h) term. Centering instead at βh,D = βh,D − (log 2)/ah gives the limiting cdf
exp{−2e−z }. This agrees with the smooth-kernel calibration of Liu and Wu (2010, Theorem 2.4) under independence and negligible bias: their λK and K2 equal R(K) and cK , and their rescaled bandwidth is h/L. Here K ⋆ (±1) = 0. Inverting the studentized limit gives the usual simultaneous interval band. For d ≥ 2, CD agrees with the smooth χd -field normalization in Qiao (2021a, Corollary 3.1). When p and Σ are constant near D, the Gaussian section has independent stationary coordinates; local whitening gives the same leading constant when they vary. Proksch (2016) constructs nonparametric bands for a scalar response with multidimensional fixed design, rather than a joint vector-norm tube. Liu et al. (2016) construct exact simultaneous confidence tubes for Gaussian multivariate linear regression with a common error covariance. Their ellipsoidal cross-sections have the same form as our Euclidean tube, but their critical value is obtained from a largest-root distribution. That finite-sample calibration is not a special case of the shrinking-bandwidth Gumbel law. The proof of Proposition 4.3 is given in the supplementary material. Returning to M, we now estimate Ω at the required relative accuracy.
21
4.3
Covariance estimation
The oracle statistic contains the unknown covariance through Ω (x) = R (K) Σ (x) /p (x). Directly averaging residual outer products with the undersmoothed inference bandwidth h can be unstable even when the numerator already has enough local observations for Gaussian approximation. We therefore use a mean-pilot bandwidth b and a covariance bandwidth g, neither of which is tied to h. The computation has three stages. We first predict the mean using observations outside each held-out fold and use the prediction errors as residuals. Next, we estimate the covariance field by smoothing their transported outer products at bandwidth g and stabilize small eigenvalues. Finally, we aggregate these covariance estimates with the squared weights of the original regression estimator at bandwidth h. For the first stage, split the sample deterministically into a fixed number K0 ≥ 2 of folds with sizes nk such that, for some cF > 0 and all large n,
min k
nk ≥ cF , n
min k
n − nk ≥ cF . n
For fold k, fit the estimator in (2.2) at bandwidth b using only observations outside that fold, with (−k)
the training-sample size in its normalization. Denote this fit by Vbb
. If observation i belongs to
fold k (i), define (−k(i))
cf ηbi,b = Yi − Vbb
(Xi ) ∈ TXi M.
(4.16)
Since observation i is not used in its own prediction, the fitted mean is independent of its held-out error conditional on the training data and Xi . For the second stage, transport residuals from nearby source points into the same target fibre and average their outer products. Smoothing at bandwidth g gives b (x) = Σ b,g
1
n X
ng d pb†g (x) i=1
n
cf Wg (x, Xi ) PXi →x ηbi,b
o⊗2
.
(4.17)
Here pb†g is the floored density estimate from Section 2, evaluated at bandwidth g. All operators in the average act on Tx M, so the estimate is defined separately at each location. If few residuals receive appreciable weight, this local matrix can be singular. We therefore floor its eigenvalues on a scale that vanishes with sample size. Choose the strictly positive intrinsic 22
residual scale (
sbcf n,b =
n 2 1 X cf ηbi,b Xi dn i=1
)
sn = n−1 ,
∨ sn ,
and set γn = {log (en)}−2 . For a positive semidefinite operator A =
j λj ej ⊗ e j
P
and a scalar
s > 0, let Fγ,s (A) =
d X
λj ∨ γ
j=1
b fl = F Set Σ b,g γn ,b scf
n,b
tr (A) +s d
ej ⊗ e j .
(4.18)
b Σ b,g . The functional calculus in (4.18) is frame-invariant. On every finite
sample its smallest eigenvalue is at least γn sn > 0 and its condition number is at most d/γn , while Proposition 4.4 shows that it is uniformly inactive with probability tending to one. For the final stage, recall that the covariance of a weighted sum uses squared weights. We evaluate the stabilized covariance field at each source point, transport it to the target, and combine it with the original regression weights. Define these weights and the resulting covariance estimate by a†i,h (x) =
Wh (x, Xi ) nhd pb†h (x)
d b Ω n,h|g (x) = nh
,
n n X †
o2
b fl (Xi ) Px→X . PXi →x Σ b,g i
(4.19)
On the event that the denominator floor is uniformly inactive, a†i,h = Wh /
j Wh , and the oracle
ai,h (x)
i=1
P
version d ΩX n,h (x) = nh
n n X †
o2
ai,h (x)
PXi →x Σ (Xi ) Px→Xi
i=1
is exactly the conditional covariance of the scaled noise term in Vbh (x). Thus (4.19) retains both the realized weights and the spatially varying error covariance. To state the accuracy needed for this replacement, let Lu = log (1/u) and write the uniform pilot-error rate as rn,b = b2 +
p
n−1 b−d log n + n−1 b−d log n. The bandwidth conditions are
b, g ↓ 0, Lh
s 2
h +
nbd → ∞, log n
Lh Lh + + g 2 + (1 + rn,b ) nhd nhd
23
ng d → ∞, Lg s
Lg Lg + ng d ng d
!
2 + rn,b → 0.
(4.20)
These conditions allow a wider neighborhood for covariance estimation than for regression, increasing the information available to estimate Σ. A concrete compatible choice, with α in (2.4), is b ≍ g ≍ (log n/n)1/(d+4) ,
1 1 <α< . d+4 d
h = n−α ,
Thus the usual smooth-function bandwidths for the nuisance fits are compatible with the undersmoothed regression bandwidth. The general conditions in (4.20) impose no fixed ratio between g and h. Proposition 4.4. Under Assumptions 2.1–2.4, (4.20), and the floor sequences specified above, (−k)
max sup Vbb
1≤k≤K0 x∈M
(x) − V (x)
x
= OP (rn,b ) .
Moreover, with n
o
b fl (x) − Σ (x) Σ (x)−1/2 eΣ,n = sup Σ (x)−1/2 Σ b,g x∈M
s
( 2
eΣ,n = OP g + (1 + rn,b )
Lg Lg + d d ng ng
!
op
,
) 2 + rn,b
.
On the event eΣ,n < 1, X b (1 − eΣ,n ) ΩX n,h (x) ⪯ Ωn,h|g (x) ⪯ (1 + eΣ,n ) Ωn,h (x)
(4.21)
simultaneously for every x. Consequently, n
o
−1/2 b sup Ω (x)−1/2 Ω n,h|g (x) − Ω (x) Ω (x)
x∈M
= OP
h2 +
s
Lh Lh + + g 2 + (1 + rn,b ) d nh nhd
op
s
Lg Lg + d d ng ng
! 2 + rn,b
.
In particular, n
o
−1/2 b Lh sup Ω (x)−1/2 Ω n,h|g (x) − Ω (x) Ω (x) x∈M
op
= oP (1) .
(4.22)
The spectral floor in (4.18) is uniformly inactive with probability tending to one. b fl at the random Xi , and The Loewner transfer in (4.21) is pathwise. Hence evaluating Σ b,g
24
reusing those observations in the h-sandwich, requires no additional covariance sample split. The factor Lh ≍ a2h in (4.22) is the exact scale required to replace the oracle inverse square root in a maximum of order ah . The proof of Proposition 4.4 is in the supplementary material. Cross-fitting makes the residual– pilot cross term conditionally centered; its remaining contribution is smaller than the covariance2 smoothing error. The squared pilot error accounts for the rn,b term. Thus broad covariance
neighborhoods can stabilize studentization while h remains small enough for inference on V . We now replace the oracle covariance by its estimate and invert the resulting Gumbel law to obtain a confidence tube. Let ω n > 0 be deterministic with ω n ↓ 0. For A = [A]−1/2 = τ
X
j λj ej ⊗ ej , write
P
(λj ∨ τ )−1/2 ej ⊗ ej .
j
The final clip makes the statistic measurable even when the h-denominator gate fails and is asymptotically inactive. The feasible statistic is
Tbn = sup x∈M
h
b Ω n,h|g (x)
i−1/2 √ ωn
n
nhd Vbh (x) − V (x)
o
. x
V | = o (1) . Indeed, the Loewner The preceding relative covariance rate implies ah |Tbn − Tn,h P
bounds make the norm perturbation at most a constant times the relative covariance error multiV = O (a ). The factor a2 in the required rate is therefore essential. plied by Tn,h P h h
For the feasible law, retain the polynomial factor in the leading excursion intensity and define the unexpanded centering (0)
λh (u) = h−d
Cvec (2d)
2
u2d−2 e−u /2 , d−1
βehvec >
√
2d − 2,
(0)
λh
βehvec = 1.
The high-branch solution is unique and exists for all sufficiently small h. Theorem 4.2 (Feasible full-norm Gumbel law). Under Assumptions 2.1–2.4 and (4.20), h
n
o
i
sup P ah Tbn − βehvec ≤ z − exp −e−z z∈R
25
−→ 0.
(4.23)
The centering in (4.23) changes the finite-h location but not the limiting law. n
ah βehvec − βhvec
o
In fact,
= o (1) . The equivalence lemma and its proof are given in the supplementary
material. To turn this limit into confidence regions, take its (1 − α) quantile:
cn,1−α = βehvec +
q1−α = − log {− log (1 − α)} ,
q1−α . ah
(4.24)
The feasible (1 − α) simultaneous confidence tube is √
Cn,1−α (x) = u ∈ Tx M :
nhd
h
b Ω n,h|g (x)
i−1/2 n ωn
u − Vbh (x)
o
≤ cn,1−α .
(4.25)
x
Theorem 4.2 implies
P {V (x) ∈ Cn,1−α (x) for every x ∈ M} −→ 1 − α.
(4.26)
At each location, (4.25) gives an ellipsoid whose orientation and relative axes come from the estimated covariance; cn,1−α supplies the common multiplicity correction across the domain. The proof of Theorem 4.2 is in the supplementary material. Algorithm 1 Simultaneous tangent-field inference at supplied bandwidths Require: Tangent observations, known geometry, kernel, level 1 − α, fixed folds, h, b, g, and vanishing numerical floors. 1: Compute the volume-corrected h-weights, density and ratio estimate (2.2). Record all unsupported targets. 2: Fit each out-of-fold source pilot at b and form residuals (4.16). Smooth transported residual outer products at g using (4.17). −2 3: Compute the residual scale with anchor n−1 and apply (4.18) with γn = {log (en)} . Insert source covariances in the realized squared-weight sandwich (4.19), then apply the final inversesquare-root clip. √ (0) 4: Compute ah and Cvec . Bracket the centering root λh (u) = 1 on u > 2d − 2. 5: if the high root does not exist or a required computation fails then 6: Report the calibration unavailable; do not substitute the low root. 7: else 8: Solve on the decreasing branch and set cn,1−α = u + q1−α /ah as in (4.24). Return (4.25). 9: end if The deterministic-bandwidth guarantee requires (2.3) and (4.20), with h < min (1, ιM /2) and b, g < ιM . For a smoothed target, the target-local covariance in Section 6 replaces the source26
(0)
covariance steps. The probability-specific equation λh (c) = − log (1 − α) on its high branch gives an asymptotically equivalent, but not identical finite-sample, cutoff (see the supplementary material). The supplementary material gives the safeguards and their coverage conditions. Numerical
. approximation of the continuum maximum preserves the limit when its error is oP a−1 h
4.4
A geometry-aware bandwidth family
The three smoothing steps need different amounts of local information, but their bandwidths can be organized around one reference radius. We describe a common choice for uniform designs on S 2 , T2 , S 3 , T3 and SO (3), allowing spatially varying covariance. Let ωd be the Euclidean unit-ball volume and set ℓM = {Vol (M) /ωd }1/d . Write κM for a bound on absolute sectional curvature. To keep the kernels inside normal neighborhoods and limit volume distortion, use √ RM = min {0.9ιM , π/ (2 κM )} ,
HM = min {0.45ιM , RM } ,
where the curvature bound contributes no cap when κM = 0. On these manifolds the normalcoordinate volume density is a radial function ΘM (t). The effective fraction of the sample used by a volume-corrected local mean at radius r is {K ⋆ (u)}2 du, ∥u∥≤1 ΘM (r∥u∥)
Z
J2,M (r) =
eM (r) =
rd . Vol (M) J2,M (r)
(4.27)
Thus neM (r) in (4.27) accounts for the unequal kernel weights when measuring the information available within radius r. Choose the smallest positive solution of neM (rn ) (rn /ℓM )4 = 1,
0 < rn ≤ RM ,
(4.28)
and use rn = RM if there is no solution in this range. The fourth power represents squared secondorder smoothing bias, whereas {neM (r)}−1 represents variance. Equation (4.28) balances these dimensionless reference scales using geometric quantities. The multiplier below allows their overall
27
scale to be adjusted. For a positive multiplier θ, set n
o
hn (θ) = min HM , θrn {log (en)}−2/(d+4) ,
bn (θ) = min {RM , θrn } , (4.29)
n
1/(d+4)
gn (θ) = min RM , θrn {log (en)}
o
.
The mean is undersmoothed for inference on V , the residual pilot uses the reference scale, and covariance smoothing uses a wider neighborhood. Dimension determines the powers, volume fixes the length scale, and curvature enters both the effective count and the geometric cap. The default is θ = 1. A finite candidate set
n
o
2j/4 : j = −4, . . . , 4
allows a factor-four range
of reference scales; out-of-fold prediction loss at bn (θ) can select among supported candidates. Here prediction loss is the mean squared tangent-norm error on held-out observations; empty neighborhoods must be excluded and local effective counts checked separately. Since J2,M (r) = R (K) + O r2 , the reference radius has order n−1/(d+4) and the geometric caps eventually become
inactive. For deterministic θ bounded above and away from zero, the logarithmic factors in (4.29) then give nhd+4 log (1/hn ) = O (1/ log n) and nhdn / (log n)5 → ∞; the pilot and covariance rates n also satisfy (4.20). The deterministic choices therefore satisfy the asymptotic inference conditions. Coverage after selecting θ by prediction loss requires a separate analysis.
5
Simulation studies
The theory provides an asymptotic simultaneous guarantee. We now examine how its analytic calibration behaves at the available sample sizes, using five geometries: S 2 and T2 in dimension two, and S 3 , T3 and SO (3) in dimension three. Spheres and tori contrast curved and flat domains, while the rotation group represents orientation data. All five experiments estimate the covariance from the observations.
5.1
Common design
The designs are uniform, the regression fields are smooth, and the errors are conditionally Gaussian. All experiments use the kernel
K ⋆ (u) = cd 1 − ∥u∥2
5
1 {∥u∥ ≤ 1} ,
28
cd =
Γ (d/2 + 6) . 120π d/2
For each n ∈ {200, 400, 600, 800, 1000, 2000, 4000}, we use 500 replicates, five cross-fitting folds and a nominal coverage of 95%. The bandwidths h, b, g were fixed after independent pilot and validation runs and before the 500 production replicates. Residuals and covariance estimates are computed by the procedure in Section 4.3. Details of the cutoff calculation and numerical maximization are given in the supplementary material. The results below use these fixed schedules; the family in Section 4.4 is a separate bandwidth construction. Table 1: Analytic-cutoff frequencies at nominal level 0.95 with fixed bandwidths, from 500 replicates per cell. Columns identify the manifold and intrinsic dimension. n d
S2 2
T2 2
S3 3
T3 3
SO (3) 3
200 400 600 800 1000 2000 4000
0.696 0.970 0.982 0.970 0.956 0.960 0.956
0.820 0.926 0.942 0.958 0.970 0.960 0.950
0.904 0.936 0.964 0.962 0.956 0.974 0.962
0.536 0.850 0.910 0.942 0.944 0.964 0.956
0.756 0.902 0.936 0.948 0.960 0.966 0.942
The analytic calibration stabilizes as the sample size increases (Table 1). For n ≥ 1000, frequencies range from 0.942 to 0.974, close to the nominal 0.95 on the reported Monte Carlo scale. At n = 200, the much lower frequencies in several geometries show the difficulty of simultaneous calibration with sparse local information. All comparisons use numerically approximated maxima. Figures 2 and 3 show the first n = 4000 replicate for the torus and sphere, respectively. Each figure restricts the fitted global tube to a geodesic, preserving its covariance estimates and global cutoff. The ellipses display how uncertainty changes along the path.
29
Figure 2: A geodesic section of the nominal 95% confidence tube on the flat torus, using the first n = 4000 production replicate. Left: the blue path γ (t) = (t, 0) with the displayed evaluation points. Right: matching blue ellipses bound the two-dimensional confidence regions in Tγ(t) T2 , translated to zero and arranged by t in the orthonormal frame e1 (t) = ∂θ1 , e2 (t) = ∂θ2 ; the red curve is Vbh {γ (t)} − V {γ (t)}. All cross-sections use one global analytic cutoff. The endpoints t = −π and t = π coincide on the torus. The embedding is schematic; estimation uses the flat product metric.
Figure 3: An equatorial section of the nominal 95% confidence tube on S 2 , using the first n = 4000 production replicate. Left: the blue equator γ (t) = (cos t, sin t, 0) and the displayed locations. Right: matching blue ellipses are zero-centered confidence cross-sections in the parallel frame e1 (t) = (− sin t, cos t, 0), e2 (t) = (0, 0, 1); the red curve is the estimation error in that frame. Every ellipse uses the same global cutoff 4.700. The endpoints t = −π and t = π represent the same fibre.
6
Randomized reconstruction of planetary-scale 500-hPa wind
The simulations assess inference for known fields. We now ask how accurately a spatially smoothed mean wind can be reconstructed from randomized space–month observations of a fixed climatological data set. The uncertainty comes from this randomized sampling, conditional on the observed wind fields. We use the NCEP–NCAR Reanalysis 1 eastward and northward monthly mean winds at 500 hPa for January 1991–December 2020, distributed by NOAA’s Physical Sciences Laboratory on 30
Figure 4: Geodesic sections of the nominal 95% analytic tube for the smoothed mean. Blue ellipses are zero-centered error sets and the red curve is the realized reconstruction error. The surface interpolates 17 evaluated fibres along each path. All rows share the physical scale and the expanded cutoff 5.692. a 2.5◦ latitude–longitude grid (Kalnay et al., 1996). Each month’s tangent field is projected by equal-area weighted least squares onto the 16 real electric and magnetic vector spherical harmonics of degrees one and two. Write Ut (x) for the resulting smooth field and
V0 (x) =
360 1 X Ut (x) , 360 t=1
Σ (x) =
360 1 X {Ut (x) − V0 (x)}⊗2 . 360 t=1
Conditional on these fields, draw n = 4000 independent uniform locations Xi ∈ S 2 and independent uniform month indices Ti ∈ {1, . . . , 360}, and observe Yi = UTi (Xi ). Then the pairs are iid, p = 1/ (4π), E (Yi | Xi = x) = V0 (x), and the error covariance is Σ (x). The responses retain units of metres per second.
31
Spatial resolution is part of the estimand. Using the same volume-corrected kernel as in the simulations, define R
Vh (x) =
h (x, q) Pq→x V0 (q) d Vol (q) S 2 WR S 2 Wh (x, q) d Vol (q)
.
The recorded choices are h = ch n−1/5 with ch = 3, 3.5, 4, corresponding to 32.72◦ , 38.17◦ , and 43.63◦ at n = 4000. The primary analysis uses ch = 3. The other choices give smoother targets by averaging over wider neighborhoods. Each resolution is fixed before calibration. We use five folds assigned by observation index modulo five and fit each training mean at the same bandwidth h as the final mean. The covariance is centered at the target after transport: P b TL (x) = Ω n,h
n
(−k(i))
2 PXi →x Yi − Vbh i Wh (x, Xi )
nh2
n
pb†h (x)
o2
(x)
o⊗2
.
(6.1)
The target-local covariance in (6.1) includes variation of the transported mean within the neighborhood as well as measurement variation. The supplementary material proves shrinking-bandwidth inference for Vh without undersmoothing. It also gives a separate fixed-resolution Gaussian result with multiplier calibration on a fixed finite grid. The reconstruction uses the expanded analytic cutoff βhvec + q0.95 /ah = 5.692375 at the primary scale. All displayed ellipsoid widths use this cutoff. It is the expanded counterpart of the unexpanded-centering cutoff in Algorithm 1. Evaluation on nested icosahedral grids of 10242 and 40962 points gives a finest-grid maximum studentized reconstruction error of 4.456, below 5.692. The fifth percentile of the effective-neighbor count is 85.8, and the covariance floor is inactive on the evaluation grid. These diagnostics describe one randomized reconstruction at the stated resolution. Figure 4 shows three geodesic sections of this single global construction. The cross-sections vary in size because the estimated covariance varies over the sphere, although all use one cutoff. Each path has length 45◦ . From top to bottom, the endpoints are the North Pole and (45◦ E, 45◦ N); (144.7◦ W, 36.4◦ N) and (136.2◦ W, 7.9◦ S); and (172.9◦ E, 23.1◦ N) and (159.6◦ E, 67.4◦ N). The outer semi-axis changes from 0.49 to 1.75 m s−1 on the first path, from 2.10 to 0.56 on the second, and from 1.43 through a maximum of 2.07 to 1.36 on the third.
32
Concluding remarks The proposed kernel estimator and confidence tube provide intrinsic estimation and simultaneous inference for tangent vector fields. Parallel transport permits local averaging across tangent spaces, while a Gaussian reference section supplies one critical value for the entire field. Estimating the covariance locally allows the tube’s orientation and width to vary across the manifold. The simulations show improving analytic calibration with increasing sample size, and the wind reconstruction illustrates the spatial variation in these uncertainty regions. The theory assumes a known compact manifold without boundary, independent observations and smooth uniformly nondegenerate covariance. Small local sample sizes can still limit the accuracy of the asymptotic calibration. For chronological wind data, a natural next step is to replace the independent-error covariance by a transported long-run covariance and establish a Gaussian approximation under temporal dependence; the resulting spatial correlation geometry must then be reanalyzed. Bias correction could permit larger regression neighborhoods, but would require uniform estimation of covariant derivatives and a new covariance calculation for the corrected process.
Supplementary material Proofs and additional technical details are provided in the supplementary material.
References Richard Arratia, Larry Goldstein, and Louis Gordon. Two moments suffice for poisson approximations: The chen–stein method. The Annals of Probability, 17(1):9–25, 1989. Hadi Beik Mohammadi, Søren Hauberg, Georgios Arvanitidis, Nadia Figueroa, Gerhard Neumann, and Leonel Rozo. Neural contractive dynamical systems. In International Conference on Learning Representations, 2024. Peter J. Bickel and Murray Rosenblatt. On some global measures of the deviations of density function estimates. The Annals of Statistics, 1(6):1071–1095, 1973.
33
Li Cai, Lijie Gu, Qihua Wang, and Suojin Wang. Simultaneous confidence bands for nonparametric regression with missing covariate data. Annals of the Institute of Statistical Mathematics, 73: 1249–1279, 2021. Victor Chernozhukov, Denis Chetverikov, and Kengo Kato. Anti-concentration and honest, adaptive confidence bands. The Annals of Statistics, 42(5):1787–1818, 2014. Victor Chernozhukov, Denis Chetverikov, and Kengo Kato. Comparison and anti-concentration bounds for maxima of gaussian random vectors. Probability Theory and Related Fields, 162(1–2): 47–70, 2015. Victor Chernozhukov, Denis Chetverikov, Kengo Kato, and Yuta Koike. Improved central limit theorem and bootstrap approximations in high dimensions. The Annals of Statistics, 50(5): 2562–2586, 2022. Swagatam Das and Vaclav Snasel. Sharp concentration bounds for bundle-valued statistics on manifolds. arXiv:2607.10592, 2026. URL https://arxiv.org/abs/2607.10592. Accepted at ICML 2026. Krzysztof Dębicki, Enkelejd Hashorva, and Peng Liu. Uniform tail approximation of homogenous functionals of gaussian fields. Advances in Applied Probability, 49(4):1037–1066, 2017. Minjie Fan, Debashis Paul, Thomas C. M. Lee, and Tomoko Matsuo. Modeling tangential vector fields on a sphere. Journal of the American Statistical Association, 113(524):1625–1636, 2018. Wolfgang Härdle and J. S. Marron. Bootstrap simultaneous error bars for nonparametric regression. The Annals of Statistics, 19(2):778–796, 1991. Guillermo Henry and Daniela Rodriguez. Robust nonparametric regression on riemannian manifolds. Journal of Nonparametric Statistics, 21(5):611–628, 2009. Michael Hutchinson, Alexander Terenin, Viacheslav Borovitskiy, So Takao, Yee Whye Teh, and Marc Peter Deisenroth. Vector-valued gaussian processes on riemannian manifolds via gauge independent projected kernels. In Advances in Neural Information Processing Systems, volume 34, pages 17160–17169, 2021.
34
Gordon J. Johnston. Probabilities of maximal deviations for nonparametric regression function estimates. Journal of Multivariate Analysis, 12(3):402–414, 1982. Eugenia Kalnay, Masao Kanamitsu, Robert Kistler, William Collins, Dennis Deaven, Lev Gandin, Mark Iredell, Suranjana Saha, Glenn White, John Woollen, Yuejian Zhu, Muthuvel Chelliah, Wesley Ebisuzaki, Wayne Higgins, John Janowiak, K. C. Mo, Chester Ropelewski, J. Wang, Ants Leetmaa, Richard W. Reynolds, Roy Jenne, and Dennis Joseph. The NCEP/NCAR 40year reanalysis project. Bulletin of the American Meteorological Society, 77(3):437–472, 1996. John M. Lee. Introduction to Riemannian Manifolds, volume 176 of Graduate Texts in Mathematics. Springer, 2 edition, 2018. Wei Liu, Yang Han, Fang Wan, Frank Bretz, and Anthony J. Hayter. Simultaneous confidence tubes in multivariate linear regression. Scandinavian Journal of Statistics, 43(3):879–885, 2016. Weidong Liu and Wei Biao Wu. Simultaneous nonparametric inference of time series. The Annals of Statistics, 38(4):2388–2421, 2010. Bruno Pelletier. Kernel density estimation on riemannian manifolds. Statistics & Probability Letters, 73(3):297–304, 2005. Bruno Pelletier. Non-parametric regression estimation on closed riemannian manifolds. Journal of Nonparametric Statistics, 18(1):57–67, 2006. Katharina Proksch. On confidence bands for multivariate nonparametric regression. Annals of the Institute of Statistical Mathematics, 68(1):209–236, 2016. Wanli Qiao. Extremes of locally stationary gaussian and chi fields on manifolds. Stochastic Processes and their Applications, 133:166–192, 2021a. Wanli Qiao. Asymptotic confidence regions for density ridges. Bernoulli, 27(2):946–975, 2021b. Wanli Qiao and Wolfgang Polonik. Extrema of rescaled locally stationary gaussian fields on manifolds. Bernoulli, 24(3):1834–1859, 2018.
35
Daniel Robert-Nicoud, Andreas Krause, and Viacheslav Borovitskiy. Intrinsic gaussian vector fields on manifolds. In Proceedings of the 27th International Conference on Artificial Intelligence and Statistics, volume 238 of Proceedings of Machine Learning Research, pages 1306–1314, 2024. Amit Singer and Hau-Tieng Wu. Vector diffusion maps and the connection laplacian. Communications on Pure and Applied Mathematics, 65(8):1067–1144, 2012. Jonathan E. Taylor and Keith J. Worsley. Random fields of multivariate test statistics, with applications to shape analysis. The Annals of Statistics, 36(1):1–27, 2008.
36