Residual Stopping and Pointwise Inference
for Nonuniform Entropic Self-Transport
Scott Schubr
Luminary Oddities LLC
October 18, 2026
Abstract
We study observable residual stopping for the half-self-Sinkhorn iteration with iid samples from
a fixed positive smooth density on the one- or two-dimensional torus. For the normalized periodic
Gaussian, assume that 𝜀 → 0, 𝑛𝜀3𝑑/2 → ∞ , and 𝑛𝜀 𝑑/2+3 → 0. Starting from the constant scaling,
we stop at the first projective square-root-mass residual below 𝑡𝑛 = 𝑛−1/2𝜀−𝛾, where 0 < 𝛾 < 𝑑 /4 is
fixed. The stopped score has the same pointwise Gaussian limit, centered at the true density score, as
the exactly scaled empirical estimator. Its attained residual lies between 𝑐𝑡𝑛 and 𝑡𝑛 with probability
tending to one for every fixed0 < 𝑐 < 1/2, although √𝑛 𝑡𝑛 → ∞. The stopping index is asymptotic to
(𝑑/4 − 𝛾) log(1/𝜀)/log 2 in probability. For the coefficient-one threshold, uniform sampling already
gives a different score limit when 𝛾 ≥ 𝑑/4. The proof retains the iteration-dependent empirical-norm
error and uses a normalized-row cancellation to control nonuniform initialization in dimension two.
Keywords: entropic self-transport; Sinkhorn iteration; statistical stopping; score estimation; pointwise
inference.
1 Introduction
Symmetric kernel scaling gives an empirical entropic self-transport score through its smooth off-sample
extension. A numerical solver is normally terminated using a normalization residual, whereas statistical
inference depends on the error in a spatial derivative. This paper identifies a family of residual tolerances
that preserves the pointwise limiting law along the actual sampled iteration.
The tolerance can exceed 𝑛−1/2 by a diverging factor. We prove this for a first-hit rule checked from
the initial scaling, rather than for a prescribed iteration number. We also determine the leading stopping
count, control the residual attained at that hit, and give a counterexample at the upper tolerance exponent.
The same observations determine the scale, iterates, residual and queried score.
1.1 Related work
Kokot’s dissertation [2], Chapter 3, gives Gaussian-RKHS scaling differentiability, the self-diagonal plus-
resolvent derivative, and a fixed-regularization Hilbert central limit theorem. Its Section 4.4.1 also gives
first-order population matching after two symmetric updates. We use these established mechanisms and
provide the joint bandwidth–sample–iteration estimates needed for residual stopping. Mordant [ 3] stud-
ies decreasing-regularization self-transport score asymptotics and the half-self iteration in a Euclidean
quadratic-cost model. Our statistical limit is proved for the exact periodic Gaussian.
Wormell and Reich [ 4] analyze empirical bistochastic normalization, resolvents and spectral conver-
gence, including an alternative symmetric scaling iteration. Cheng and Landa [ 5] show that approximate
nonuniform normalization at statistical accuracy suffices for graph-Laplacian consistency. Maeda, Y ao
and Nitanda [7] combine sampling and iteration errors for two-sample Schrödinger bridges. Spectral dis-
crepancy stopping is developed in, among other works, Celisse and Wahl [ 6]. These results establish the
general statistical motivation for early termination.
1

The present conclusions concern the specific nonlinear first residual hit: its admissible exponent range,
attained residual and leading count, together with preservation of the exact estimator’s true-score limit.
The earlier residual-certificate paper [ 8] gives bounds for arbitrary positive candidates and empirical
sharpness witnesses. It does not give the actual-trajectory stopping theorem proved here. The bandwidth
windows below are sufficient; no optimal sampling range or wall-clock complexity claim is made.
2 Model and main theorem
Let 𝑑 ∈ { 1, 2} and let 𝑓 ∈ 𝐶6(T𝑑) be a fixed probability density with min 𝑓 > 0. Set 𝜇 = 𝑓 𝑑𝑥 and let
𝜇𝑛 = 𝑛−1 Í𝑛
𝑖=1 𝛿𝑋𝑖 for iid 𝑋𝑖 ∼ 𝜇. Lebesgue measure on the unit torus has total mass one. The kernel is
𝑘 𝜀 (𝑧) =
Õ
ℓ ∈Z𝑑
(2𝜋𝜀) −𝑑/2𝑒− | 𝑧+ℓ | 2/( 2𝜀) , 𝜆 𝑘 = 𝑒−2 𝜋2 𝜀 | 𝑘 | 2
. (1)
We also write 𝑘 𝜀 (𝑥, 𝑦 ) = 𝑘 𝜀 (𝑥 − 𝑦). Here 𝜆𝑘 are the Fourier multipliers and 𝐻 𝜀 = exp(𝜀Δ/2) is convo-
lution against Lebesgue measure. Write 𝐾𝜈 𝑏(𝑥) =
¯
𝑘 𝜀 (𝑥 − 𝑦)𝑏(𝑦) 𝑑𝜈 (𝑦) and 𝐾𝑛 = 𝐾𝜇𝑛. Throughout,
𝜀 = 𝜀𝑛 → 0, 𝑛𝜀 3𝑑/2 → ∞, 𝑛𝜀 𝑑/2+3 → 0. (2)
Constants may depend on this fixed 𝑓 and on fixed exponents. There is no assertion uniform over an
unbounded class of density derivative norms.
Let b𝑎 be the positive empirical self-scale, b𝑎𝐾 𝑛b𝑎 = 1. At the sample sites it exists uniquely: the strictly
convex coercive function
Φ(𝑢) = 1
2𝑛2
Õ
𝑖,𝑙
𝑘 𝜀 (𝑋𝑖 − 𝑋𝑙)𝑒𝑢𝑖+𝑢𝑙 − 1
𝑛
Õ
𝑖
𝑢𝑖
has a unique minimizer and its gradient equations are the scaling equations. Coercivity and strict convex-
ity follow already from the positive diagonal terms. Extend the scale by b𝑎(𝑥) = 1/𝐾𝑛b𝑎(𝑥) away from the
sites. Starting with 𝑏0 = 1, use the averaged symmetric iteration [ 1, 3]
𝑏 𝑗+1 =

𝑏 𝑗/( 𝐾𝑛𝑏 𝑗)
	1/2 , 𝑠 𝑗 = 2∇ log 𝐾𝑛𝑏 𝑗, b𝑠 = 2∇ log 𝐾𝑛b𝑎. (3)
Only site values are needed for each convolution. Define
𝐹𝑗 = (𝑏 𝑗 𝐾𝑛𝑏 𝑗)1/2, 𝑆 𝑗 = 𝜇𝑛𝐹𝑗, 𝑟 𝑗 =


𝐹𝑗/𝑆 𝑗 − 1



𝑛 , ∥𝑢∥2
𝑛 = 𝜇𝑛 | 𝑢| 2 . (4)
Both 𝑟 𝑗 and 𝑠 𝑗 are invariant under positive scalar rescaling of 𝑏 𝑗. For fixed 𝛾 > 0, put
𝑡𝑛 = 𝑛−1/2𝜀−𝛾, 𝐿 𝑛 = max{ 3, ⌈4 log 𝑛⌉} ,
𝜏𝛾 = min { 𝑗 ∈ { 0, . . . , 𝐿𝑛} : 𝑟 𝑗 ≤ 𝑡𝑛} ∪ { 𝐿𝑛} . (5)
The budget makes the rule finite on every sample. Set
𝐴𝑛 = √𝑛 𝜀𝑑/4+1/2, 𝐽 𝑑 =
¹
R𝑑
(2𝜋𝜉1)2
(1 + 𝑒2 𝜋2 | 𝜉 | 2
)2
𝑑𝜉 > 0. (6)
Theorem 2.1. Assume (2) and fix 0 < 𝛾 < 𝑑 /4. Then
𝜏𝛾
log(1/𝜀) − →𝑑/4 − 𝛾
log 2 in probability. (7)
For every fixed 0 < 𝑐 < 1/2,
P{ 𝜏𝛾 < 𝐿 𝑛, 𝑐𝑡 𝑛 ≤ 𝑟 𝜏𝛾 ≤ 𝑡𝑛} → 1. (8)
2

For each fixed 𝑥 ∈ T𝑑,
𝐴𝑛


𝑠𝜏𝛾 (𝑥) − b𝑠(𝑥)


 → 0 in probability, (9)
𝐴𝑛{ 𝑠𝜏𝛾 (𝑥) − ∇ log 𝑓 (𝑥)} ⇒ 𝑁𝑑

0, 4𝐽𝑑
𝑓 (𝑥) 𝐼𝑑

. (10)
The exact empirical score b𝑠 satisfies the same centered limit.
In particular √𝑛 𝑟𝜏𝛾 → ∞ in probability. For 𝜀 = 𝑛−𝜅 the sufficient power windows are
𝑑 = 1 : 2 /7 < 𝜅 < 2/3, 𝑑 = 2 : 1 /4 < 𝜅 < 1/3.
For fixed 0 < 𝛾 1 < 𝛾 2 < 𝑑 /4, subtraction on the same sample gives
𝜏𝛾1 − 𝜏𝛾2
log(1/𝜀) − →𝛾2 − 𝛾1
log 2 in probability. (11)
This compares update counts, with no claim about cost per update. The upper exponent obstruction is
proved in Section 7.
3 Population geometry and the empirical scaling branch
Throughout this section, 𝑑 ∈ { 1, 2} , 𝑓 is the fixed density in the main theorem, and 𝜇 = 𝑓 𝑑𝑥 . Constants
can depend on 𝑓 and 𝑑. Write 𝐻 𝜀 for convolution against the kernel in ( 1) with respect to Lebesgue
measure and write
𝐾𝜇𝑢 = 𝐻 𝜀 ( 𝑓 𝑢), 𝐾 𝑛𝑢 = 1
𝑛
𝑛Õ
𝑖=1
𝑘 𝜀 ( · , 𝑋𝑖)𝑢(𝑋𝑖).
Let 𝑎 denote the positive population scale, 𝑎𝐾 𝜇𝑎 = 1, and put
𝑣 = 𝑎−1 = 𝐾𝜇𝑎, 𝑤 = 𝐻 𝜀 𝑓 , 𝑞 = 𝑎
p
𝑓 , 𝑑 1 = log 𝑣 − 1
2 log 𝑤.
All bounds below hold as 𝜀 ↓ 0.
Lemma 3.1 (Population approximation). The positive population scale is unique and satisfies
∥𝑞 − 1∥𝐶2 ≤ 𝐶𝜀, ∥𝑑1∥𝐶2 ≤ 𝐶𝜀, ∥ − 2∇ log 𝑎 − ∇ log 𝑓 ∥∞ ≤ 𝐶𝜀. (12)
In particular, 𝑎, 𝑣, 𝑤 are bounded above and below, and their 𝐶2 norms are bounded uniformly for small
𝜀.
Proof. Set 𝑔 𝑓 =
p
𝑓 . The scaling equation becomes
F𝜀 (𝑞) := 𝑞𝐻 𝜀 (𝑞𝑔 𝑓 ) − 𝑔 𝑓 = 0.
We solve it in the Banach algebra 𝐶2(T𝑑). First, ∥𝐻 𝜀𝑔 𝑓 − 𝑔 𝑓 ∥𝐶2 ≤ 𝐶𝜀 ∥𝑔 𝑓 ∥𝐶4 by the heat semigroup
identity. Second, (𝐼 + 𝐻 𝜀) −1 is uniformly bounded on 𝐶2. Indeed, its Fourier multiplier is 1 − 𝐹 (√𝜀𝑘 ),
where
𝐹 (𝜉) = 1
1 + exp(2𝜋2| 𝜉 | 2) .
The function 𝐹 is Schwartz. Periodization of its rescaled inverse Fourier transform gives a convolution
kernel whose 𝐿1 norm is bounded independently of 𝜀. The same bound holds after commuting derivatives
through convolution.
At 𝑞 = 1,
𝐷F𝜀 (1)𝑢 = (𝐻 𝜀𝑔 𝑓 )𝑢 + 𝐻 𝜀 (𝑔 𝑓 𝑢) = 𝑔 𝑓 (𝐼 + 𝐻 𝜀)𝑢 + 𝐸 𝜀𝑢,
3

where ∥𝐸 𝜀 ∥𝐶2→𝐶2 = 𝑂 (√𝜀). To check the last bound, differentiate the commutator [𝐻 𝜀, 𝑀𝑔 𝑓 ]𝑢 through
order two. Each resulting term is a commutator with a derivative of 𝑔 𝑓 ; its kernel contains the difference
of that derivative at two points. The Gaussian first absolute moment is 𝑂 (√𝜀), so each term is bounded
by 𝐶√𝜀∥𝑢∥𝐶2. The multiplication error 𝐻 𝜀𝑔 𝑓 − 𝑔 𝑓 has the stronger 𝑂 (𝜀) bound. Consequently 𝐷F𝜀 (1)
has a uniformly bounded inverse for small 𝜀. The remaining term in F𝜀 (1 + 𝑢) is 𝑢𝐻 𝜀 (𝑔 𝑓 𝑢), bounded by
𝐶 ∥𝑢∥2
𝐶2 with the corresponding local Lipschitz estimate. A contraction on a ball of radius 𝐶𝜀 therefore
yields a positive solution with 𝑞 = 1 + 𝑂𝐶2 (𝜀).
For uniqueness, two positive scales 𝑎, 𝑏 obey
0 ≤ ⟨ 𝑎 − 𝑏, 𝐾 𝜇 (𝑎 − 𝑏)⟩ 𝐿2 ( 𝜇) = −
¹ (𝑎 − 𝑏)2
𝑎𝑏 𝑑𝜇 ≤ 0.
Here positivity follows from the nonnegative Fourier multipliers of the kernel. Hence 𝑎 = 𝑏. Finally
𝑎 = 𝑞/
p
𝑓 , 𝑤 = 𝑓 + 𝑂𝐶2 (𝜀), and the elementary product, reciprocal, and logarithm estimates in 𝐶2 give
every conclusion in ( 12). □
Let H = H𝜀 be the real reproducing-kernel Hilbert space of 𝑘 𝜀, with feature 𝑘 𝑦 = 𝑘 𝜀 ( · , 𝑦). Define
𝑆 : 𝐿2(𝜇) − → H, 𝑆𝑢 =
¹
𝑎(𝑦)𝑘 𝑦𝑢(𝑦) 𝑑𝜇(𝑦), 𝑆∗ℎ(𝑦) = 𝑎(𝑦)ℎ(𝑦),
𝐶 = 𝑆𝑆 ∗, 𝑃 = 𝑆∗𝑆, 𝐵 = 1
2 (𝐼 − 𝐶), 𝑅 = (𝐼 + 𝐶) −1.
For a fixed query 𝑥 and component 𝜈 ∈ { 1, . . . , 𝑑} , set
𝜃 𝑥,𝜈 = 𝜕𝑥𝜈 { 𝑎(𝑥)𝑘 𝑥 } , 𝛼 𝑑 = 𝑑
4 + 1
2 .
For a vector section 𝜃 𝑥, inner products with a Hilbert-space element are taken componentwise.
Lemma 3.2 (Feature geometry). The operator 𝑃 is a self-adjoint positive Markov operator. Moreover,
∥𝑆∥ = ∥𝐶 ∥ = 1, ∥𝐵∥ ≤ 1
2 , ∥𝑅∥ ≤ 1, ∥𝑣∥ H = 1, 𝐶𝑣 = 𝑣, 𝑆 ∗𝑣 = 1, (13)
sup
𝑦
∥𝑘 𝑦 ∥ H ≤ 𝐶𝜀 −𝑑/4, sup
𝑥
∥𝜃 𝑥,𝜈 ∥ H ≤ 𝐶𝜀 − 𝛼𝑑 , tr 𝐶 ≤ 𝐶𝜀 −𝑑/2. (14)
For every fixed 𝑥 and 𝜈, ℓ ∈ { 1, . . . , 𝑑} ,
𝜀𝑑/2+1⟨𝜃 𝑥,𝜈 , 𝑅𝐶 𝑅𝜃𝑥,ℓ ⟩H − → 𝐽𝑑
𝑓 (𝑥) 1{ 𝜈=ℓ } , 𝐽 𝑑 =
¹
R𝑑
(2𝜋𝜉1)2
[1 + exp(2𝜋2| 𝜉 | 2)] 2 𝑑𝜉. (15)
Every fixed interval of positive length strictly inside (0, 1) contains at least 𝑐𝜀 −𝑑/2 eigenvalues of 𝐶,
counted with multiplicity, for all suﬀiciently small 𝜀.
Proof. The kernel of 𝑃 relative to 𝜇 is 𝑎(𝑥)𝑘 𝜀 (𝑥, 𝑦 )𝑎(𝑦). It is symmetric, positive semidefinite, and has
row sum one by the scaling equation. Jensen’s inequality gives ∥𝑃∥ ≤ 1, and 𝑃1 = 1 gives equality. The
singular-value identities for 𝑆 yield the asserted norms. Also
∥𝑣∥2
H =
¹
𝑎(𝑦)𝑣(𝑦) 𝑑𝜇(𝑦) = 1, 𝐶𝑣 =
¹
𝑎(𝑦)2𝑘 𝑦𝑣(𝑦) 𝑑𝜇(𝑦) = 𝑣,
and 𝑆∗𝑣 = 𝑎𝑣 = 1. Gaussian Fourier sums give
𝑘 𝜀 (𝑦, 𝑦) =
Õ
𝑘 ∈Z𝑑
𝑒−2 𝜋2 𝜀 | 𝑘 | 2
= 𝑂 (𝜀−𝑑/2), ∥𝜕𝑥𝜈 𝑘 𝑥 ∥2
H =
Õ
𝑘
(2𝜋𝑘 𝜈)2𝑒−2 𝜋2 𝜀 | 𝑘 | 2
= 𝑂 (𝜀−𝑑/2−1).
Together with Lemma 3.1, these prove the section bounds. The trace formula is
tr 𝐶 = 𝑘 𝜀 (0, 0)
¹
𝑎2 𝑑𝜇 = 𝑘 𝜀 (0, 0)
¹
𝑞2 𝑑𝑥 = 𝑂 (𝜀−𝑑/2).
4

There is a convenient comparison within the same Hilbert space. Put
𝐶𝐻 =
¹
𝑘 𝑦 ⊗ 𝑘 𝑦 𝑑𝑦, 𝑅 𝐻 = (𝐼 + 𝐶𝐻) −1.
Since 𝑎2 𝑓 = 𝑞2, comparison of quadratic forms gives
∥𝐶 − 𝐶𝐻 ∥op ≤ ∥ 𝑞2 − 1∥∞∥𝐶𝐻 ∥op = 𝑂 (𝜀). (16)
Here 𝐶𝐻 has eigenvalues 𝜆𝑘 = exp(−2𝜋2𝜀| 𝑘 | 2): in the complexification of H , its orthonormal Fourier
basis consists of √𝜆𝑘 exp(2𝜋𝑖𝑘 · 𝑥). In particular ∥𝐶𝐻 ∥ = 1. Resolvent identities give
∥𝑅𝐶 𝑅 − 𝑅𝐻𝐶𝐻 𝑅𝐻 ∥op = 𝑂 (𝜀).
Also 𝜃 𝑥,𝜈 = 𝑎(𝑥)𝜕𝑥𝜈 𝑘 𝑥 +( 𝜕𝑥𝜈 𝑎(𝑥)) 𝑘 𝑥, where the second term has norm𝑂 (√𝜀 𝜀 − 𝛼𝑑 ). After multiplication
by 𝜀2𝛼𝑑 , the operator replacement therefore has error 𝑂 (𝜀) and these section replacements have error
𝑂 (√𝜀). Since 𝑎(𝑥)2 → 1/ 𝑓 (𝑥), the Gram limit reduces to
𝜀𝑑/2+1
Õ
𝑘 ∈Z𝑑
(2𝜋𝑘 𝜈) (2𝜋𝑘 ℓ)
 
𝑒−2 𝜋2 𝜀 | 𝑘 | 2
1 + 𝑒−2 𝜋2 𝜀 | 𝑘 | 2
! 2
.
This is a Gaussian-dominated Riemann sum. Its limit is 𝐽𝑑 for equal components and zero for unequal
components, proving (15).
Finally, by (16) and the compact self-adjoint min–max principle, the ordered eigenvalues of𝐶 and 𝐶𝐻
differ by at most 𝐶𝜀 . Inside a fixed interval 𝐼 ⊂ ( 0, 1), choose a smaller interval at positive distance from
its endpoints. The heat eigenvalues in that smaller interval correspond to lattice points in a fixed annulus
after multiplication by √𝜀; their number is bounded below by 𝑐𝜀 −𝑑/2. All their perturbed eigenvalues
belong to 𝐼 for small 𝜀, proving the count. □
For the empirical statements, all random quantities use the same sample. Define
𝑔 = (𝜇𝑛 − 𝜇) (𝑎𝑘 ), 𝑔 0 = (𝜇𝑛 − 𝜇) (𝑘),
𝜁 = 𝑛−1/2𝜀−𝑑/4, 𝜂 = 𝑛−1/2𝜀−𝑑/2, 𝐶 0,𝑛 =
¹
𝑎(𝑦)2𝑘 𝑦 ⊗ 𝑘 𝑦 𝑑𝜇𝑛 (𝑦).
The notation 𝑂P, H below means a bound in Hilbert-space norm.
Lemma 3.3 (Quantitative empirical branch). Under (2), the exact empirical scale b𝑎, extended by b𝑎𝐾 𝑛b𝑎 =
1, obeys
b𝑣 − 𝑣 = 𝑅𝑔 + 𝑂P, H (𝜂𝜁 ), ∥b𝑣 − 𝑣∥ H = 𝑂P(𝜁), b𝑣 = b𝑎−1, (17)
∥b𝑎 − 𝑎∥∞ + √𝜀∥∇(b𝑎 − 𝑎) ∥∞ = 𝑂P(𝜂). (18)
In addition,
∥𝑔∥ H + ∥ 𝑔0∥ H = 𝑂P(𝜁), Cov(𝑔) = 𝐶 − 𝑣 ⊗ 𝑣
𝑛 , ⟨𝑣, 𝑔⟩H = 0, (19)
∥𝐶0,𝑛 − 𝐶 ∥op ≤ ∥ 𝐶0,𝑛 − 𝐶 ∥HS = 𝑂P(𝜂). (20)
Define 𝑆𝑛, 𝑆∗
𝑛, b𝐶, 𝑃 𝑛 by replacing 𝑎, 𝜇 in 𝑆, 𝑆∗, 𝐶, 𝑃 with b𝑎, 𝜇 𝑛. Then
∥𝑆𝑛 ∥ = ∥ b𝐶 ∥ = 1, ∥ b𝐶 − 𝐶 ∥op = 𝑂P(𝜂), (21)
sup
𝑥
∥b𝜃 𝑥,𝜈 − 𝜃 𝑥,𝜈 ∥ H = 𝑂P(𝜂𝜀 − 𝛼𝑑 ), b𝜃 𝑥,𝜈 = 𝜕𝑥𝜈 { b𝑎(𝑥)𝑘 𝑥 } . (22)
On a common event of probability tending to one, evaluation satisfies ∥ℎ∥ 𝐿2 ( 𝜇𝑛 ) ≤ 𝐶 ∥ℎ∥ H for every
ℎ ∈ H .
5

Proof. The expectation of 𝑎(𝑋)𝑘 𝑋 is 𝑣, and its second moment operator is 𝐶. Hence the covariance
identity follows. Moreover ⟨𝑣, 𝑎 (𝑦)𝑘 𝑦⟩ = 𝑎(𝑦)𝑣(𝑦) = 1, so the claimed orthogonality is exact. The
feature norm bound gives E∥𝑔∥2
H + E∥𝑔0∥2
H ≤ 𝐶𝜁 2. Similarly,
E∥𝐶0,𝑛 − 𝐶 ∥2
HS ≤ 1
𝑛 E∥𝑎(𝑋)2𝑘 𝑋 ⊗ 𝑘 𝑋 ∥2
HS ≤ 𝐶𝑛 −1𝜀−𝑑 = 𝐶𝜂 2.
These are Hilbert-space second-moment bounds; no finite-dimensional operator inequality is needed.
Consider G𝑛 (𝑧) = 𝑧 − 𝐾𝑛 (1/𝑧) in a neighborhood of 𝑣 in H . The function 𝑣 is bounded below, and
G𝑛 (𝑣) = −𝑔, 𝐷 G𝑛 (𝑣) = 𝐼 + 𝐶0,𝑛, ∥ (𝐼 + 𝐶0,𝑛) −1∥ ≤ 1.
For ℎ small enough to keep 𝑣 + ℎ positive, the reciprocal remainder is
𝐾𝑛
 ℎ2
𝑣2(𝑣 + ℎ)

.
Its H norm is bounded by 𝐶𝜀 −𝑑/4∥ℎ∥2
H on an event where ∥𝐶0,𝑛 ∥ ≤ 2. Indeed, write the remainder as
the feature operator with weights 𝑎 applied to (𝑎ℎ)ℎ/(𝑣 + ℎ); its operator norm is bounded by ∥𝐶0,𝑛 ∥1/2,
while ∥𝑎ℎ ∥ 𝐿2 ( 𝜇𝑛 ) ≤ ∥ 𝐶0,𝑛 ∥1/2∥ℎ∥ H and ∥ℎ∥∞ ≤ 𝐶𝜀 −𝑑/4∥ℎ∥ H. The same argument gives the local
Lipschitz remainder bound. Since 𝜀−𝑑/4𝜁 = 𝜂 → 0, a contraction on a ball of radius 𝑂P(𝜁) yields a
positive solution and
b𝑣 − 𝑣 = (𝐼 + 𝐶0,𝑛) −1𝑔 + 𝑂P, H (𝜂𝜁 ).
For precision, one can first restrict to ∥𝑔∥ ≤ 𝑀 𝜁 for fixed 𝑀, perform the contraction for all sufficiently
large 𝑛, and then let 𝑀 increase. The resolvent identity and (20) replace (𝐼 +𝐶0,𝑛) −1 by 𝑅 at cost 𝑂P, H (𝜂𝜁 ).
The positive solution is the exact empirical branch: at sample sites its reciprocal satisfies the finite
scaling equations, whose positive solution is unique by the same positive-semidefinite energy identity as
in Lemma 3.1. Off sample, b𝑣 = 𝐾𝑛b𝑎 fixes the stated extension. There is no extension to arbitrary signed
measures in this argument.
The feature and derivative-feature bounds applied to ℎ = b𝑣 − 𝑣 show
∥ℎ∥∞ + √𝜀∥∇ ℎ∥∞ = 𝑂P(𝜂).
Taking reciprocals, using the deterministic lower bound on 𝑣, proves (18). The empirical operator 𝑃𝑛 is
positive, self-adjoint, and Markov on 𝐿2(𝜇𝑛); hence the norm identities in (21) follow as in the population
case. Also b𝑎/𝑎 = 1 + 𝑂P(𝜂) uniformly. Comparing quadratic forms gives
∥ b𝐶 − 𝐶0,𝑛 ∥op ≤ 𝑂P(𝜂) ∥𝐶0,𝑛 ∥op,
which proves the remaining operator bound. Differentiating (b𝑎 − 𝑎) (𝑥)𝑘 𝑥 and using ( 18) proves ( 22).
Finally, since 𝑎 is bounded below,
∥ℎ∥2
𝐿2 ( 𝜇𝑛 ) ≤ 𝐶
¹
𝑎2ℎ2 𝑑𝜇𝑛 = 𝐶 ⟨ℎ, 𝐶0,𝑛 ℎ⟩H ≤ 𝐶 ∥ℎ∥2
H
on the same event. This proves the simultaneous evaluation assertion. □
4 Initialization and the sampled trajectory
All operators below use the same observations. Write 𝑒 𝑗 = log(𝑏 𝑗/b𝑎) and 𝐵𝑛 = (𝐼 −𝑃𝑛)/2, with empirical
operators as in Lemma 3.3. Off sample, 𝑃𝑛𝑢(𝑥) = b𝑎(𝑥)𝐾𝑛 (b𝑎𝑢) (𝑥) is a Markov row. The exact log update
is
Ψ𝑛 (𝑒) = 1
2 { 𝑒 − log(𝑃𝑛 (exp 𝑒))} , 𝑒 𝑗+1 = Ψ 𝑛 (𝑒 𝑗). (23)
Its uniform norm is nonincreasing. For small ∥𝑒∥∞, Taylor expansion, Jensen’s inequality and the PSD
Markov property give
∥Ψ𝑛 (𝑒) − 𝐵𝑛𝑒∥𝑛 ≤ 𝐶 ∥𝑒∥∞ ∥𝑒∥𝑛 , ∥𝐵𝑛 ∥op ≤ 1
2 . (24)
6

Lemma 4.1 (T wo-update initialization). Let 𝑣 = 1/𝑎, 𝑤 = 𝐻 𝜀 𝑓 , and 𝑑1 = log 𝑣 − 1
2 log 𝑤. Put 𝜁 =
𝑛−1/2𝜀−𝑑/4, 𝜂 = 𝑛−1/2𝜀−𝑑/2 and 𝑢 = 𝑒1 − 𝑑1. Then
𝑢 = 𝑎𝑅𝑔 − 1
2 𝑔0/𝑤 + 𝑂 𝐿2 ( 𝜇𝑛 ),P(𝜂𝜁 ), (25)
∥𝑢∥𝑛 = 𝑂P(𝜁), ∥𝑢∥∞ = 𝑂P(𝜂),
∥𝐵𝑛𝑑1∥𝑛 = 𝑂P(𝜀2 + 𝜀3/2𝜂), (26)
𝑒2 = 𝐵𝑛𝑢 + 𝐸2, ∥𝐸2∥𝑛 = 𝑂P(𝜀2 + 𝜀3/2𝜂 + 𝜀𝜁 + 𝜂𝜁 ), (27)
∥𝑒2∥𝑛 = 𝑂P(𝜀2 + 𝜁), ∥𝑒2∥∞ = 𝑂P(𝜀2 + 𝜂). (28)
Proof. The first iterate is𝑏1 = (𝐾𝑛1) −1/2. Expand 𝑒1 = logb𝑣− 1
2 log 𝐾𝑛1, using Lemma 3.3, 𝐾𝑛1 = 𝑤+𝑔0,
and the uniform value errors 𝑂P(𝜂). A quadratic logarithmic remainder has sample 𝐿2 norm at most its
uniform error times its 𝐿2 error. This proves ( 25).
For (26), Lemma 3.1 gives ∥𝑑1∥𝐶2 ≤ 𝐶𝜀 . Thus
| 𝑑1(𝑦) − 𝑑1(𝑥)| ≤ 𝐶𝜀 dist(𝑥, 𝑦 ), ∥(𝐼 − 𝑃)𝑑1∥∞ ≤ 𝐶𝜀 2.
The second estimate follows from the heat-semigroup identity applied to 𝑎 𝑓 𝑑1 − 𝑑1(𝑥)𝑎 𝑓 , using bounded
𝐶2 norms of 𝑎 𝑓 . Replace the normalized row weights b𝑎(𝑦) by 𝑎(𝑦). Their relative change is 𝑂P(𝜂), so
its cost is this factor times the empirical row mean of | 𝑑1(𝑦) − 𝑑1(𝑥)| . The Gaussian-image bound
dist(𝑥, 𝑦 )𝑘 𝜀 (𝑥 − 𝑦) ≤ 𝐶√𝜀 𝑘 2𝜀 (𝑥 − 𝑦)
and the RKHS value estimate at bandwidth 2𝜀 make this mean 𝑂P(𝜀3/2), uniformly in 𝑥. Denominators
are bounded below since 𝐾𝑛𝑎 = 𝐾𝜇𝑎 + 𝑂P(𝜂) uniformly. The weight cost is therefore 𝑂P(𝜀3/2𝜂).
For the remaining fixed-𝑎 empirical row, its centered numerator is
𝑁𝑛 (𝑥) = (𝜇𝑛 − 𝜇)𝑦 { 𝑘 𝜀 (𝑥 − 𝑦)𝑎(𝑦) [𝑑1(𝑦) − 𝑑1(𝑥)]} .
At 𝑥 = 𝑋𝑖 its own diagonal vanishes. Conditional on 𝑋𝑖 = 𝑥, its mean is 𝑂 (𝜀2/𝑛) and its variance is at
most
𝐶𝜀 2
𝑛
¹
𝑘 𝜀 (𝑥 − 𝑦)2 dist(𝑥, 𝑦 )2 𝑑𝜇(𝑦) ≤ 𝐶
𝑛 𝜀3−𝑑/2.
Averaging conditional second moments gives∥𝑁𝑛 ∥𝑛 = 𝑂P(𝜀3/2𝜁 +𝜀2/𝑛). The uniform denominator error
is 𝑂P(𝜂) and multiplies a population numerator 𝑂 (𝜀2). This proves ( 26) without an independent-query
assumption.
Taylor-expand (23) at deterministic 𝑑1. The tilted row 𝑄𝑛,𝑑 1 ℎ = 𝑃𝑛 (𝑒𝑑1 ℎ)/( 𝑃𝑛𝑒𝑑1) differs from 𝑃𝑛
by 𝑂 (𝜀) in sample 𝐿2 operator norm, by bounded multiplication and the contraction of 𝑃𝑛. Its quadratic
remainder in 𝑢 has norm at most 𝐶 ∥𝑢∥∞ ∥𝑢∥𝑛. Also ∥Ψ𝑛 (𝑑1) − 𝐵𝑛𝑑1∥𝑛 ≤ 𝐶𝜀 2. These give ( 27) and
its 𝐿2 consequence. In the uniform norm Ψ𝑛 is Lipschitz with constant one. The coarser row estimate
∥Ψ𝑛 (𝑑1)∥∞ = 𝑂P(𝜀2 + 𝜀𝜂) uses ∥𝑑1∥∞ 𝑂P(𝜂) for weight and empirical-row errors. With ∥𝑢∥∞ = 𝑂P(𝜂)
this proves (28). □
Lemma 4.2 (Common trajectory and score proxies). There are nonnegative𝐸𝑛 = 𝑜P(1) and 𝑀𝑛 = 𝑂P(1),
common to all integers 𝑗 ≥ 2, such that
𝑒 𝑗 = 𝑆∗
𝑛𝐵 𝑗 𝑅𝑔 + 𝐸 𝑗,𝑛 , √𝑛


𝐸 𝑗,𝑛



𝑛 ≤ 𝐸𝑛 + 𝑀𝑛𝜀1/2−𝑑/42− 𝑗, (29)
𝐴𝑛


𝑠 𝑗 (𝑥) − b𝑠(𝑥) − 𝑉 𝑗 (𝑥)


 ≤ 𝐸𝑛 + 𝑀𝑛𝜀1/2−𝑑/42− 𝑗, 𝑉 𝑗 (𝑥) = 2

𝜃 𝑥, 𝐶 𝐵𝑗 𝑅𝑔

. (30)
The score remainder is uniform in 𝑥; vector-section inner products are componentwise. At each fixed 𝑥,
for every deterministic 𝐽𝑛 → ∞,
sup
𝑗 ≥ 𝐽𝑛
𝐴𝑛


𝑠 𝑗 (𝑥) − b𝑠(𝑥)


 → 0 in probability. (31)
7

Proof. Let 𝛿𝑛 = ∥𝑒2∥∞. On one event of probability tending to one, (24) and uniform-norm nonexpansion
imply


𝑒 𝑗



𝑛 ≤ ( 3/4) 𝑗 −2 ∥𝑒2∥𝑛 for all 𝑗 ≥ 2. Iterating the remainder yields
sup
𝑗 ≥2



𝑒 𝑗 − 𝐵 𝑗 −2
𝑛 𝑒2




𝑛
= 𝑂P{( 𝜀2 + 𝜂) (𝜀2 + 𝜁)} . (32)
Every error in ( 27) and ( 32) vanishes after multiplication by √𝑛 under (2):
√𝑛𝜀2 → 0, √𝑛𝜀3/2𝜂 = 𝜀3/2−𝑑/2 → 0, √𝑛𝜂𝜁 = (𝑛𝜀3𝑑/2) −1/2 → 0.
The remaining terms are bounded by 𝜀1−𝑑/4, 𝜀2−𝑑/4, 𝜀2−𝑑/2 and √𝑛𝜀4, all tending to zero.
The commutator
𝑞𝑛 (𝑥) = 𝑔0(𝑥) − 𝑣(𝑥)𝑔(𝑥) = (𝜇𝑛 − 𝜇)𝑦 𝑘 𝜀 (𝑥 − 𝑦) [1 − 𝑣(𝑥)𝑎(𝑦)]
has zero diagonal at 𝑥 = 𝑋𝑖. Its coefficient is bounded by 𝐶 dist(𝑥, 𝑦 ), its conditional bias is 𝑂 (𝜀/𝑛), and
its conditional variance is 𝑂 (𝜀1−𝑑/2/𝑛). Consequently ∥𝑞𝑛 ∥𝑛 = 𝑂P(√𝜀 𝜁 ). Using 𝑣/𝑤 − 𝑎 = 𝑂 (𝜀) and
𝑅 − 1
2 𝐼 = 𝐵𝑅 in (25) gives
𝑢 = 𝑆∗
𝑛𝐵𝑅𝑔 + 𝑂 𝐿2 ( 𝜇𝑛 ),P(𝜀𝜁 + 𝜂𝜁 + √𝜀 𝜁 ).
The intertwining 𝐵𝑛𝑆∗
𝑛 = 𝑆∗
𝑛b𝐵, with b𝐵 = (𝐼 − b𝐶)/2, is exact. Both half-filters have norm at most 1/2.
Telescoping their powers contributes at most 𝐶 𝑗2− 𝑗𝜂 ∥𝑔∥ H = 𝑂P(𝜂𝜁 ) uniformly in 𝑗. The commuta-
tor retains its factor 2− 𝑗. This proves ( 29) with common random bounds. In 𝑑 = 2 its fixed- 𝑗 root-𝑛
commutator term is not negligible.
The exact score difference is 2∇ log(𝑃𝑛 (exp 𝑒 𝑗)). The derivative row is 𝑆∗
𝑛b𝜃 𝑥, of norm at most
𝑂P(𝜀−𝑑/4−1/2). Expanding the exponential and the logarithm bounds its nonlinear score remainder
by 𝐶𝜀 −𝑑/4−1/2𝛿𝑛


𝑒 𝑗



𝑛. After multiplication by 𝐴𝑛, this is bounded by the root- 𝑛 version of ( 32).
The linear proxy is 2
D
b𝜃 𝑥, b𝐶 𝐵𝑗 𝑅𝑔
E
up to the trajectory error. Replacing sections and operators costs
𝑂P(𝜀−𝑑/4−1/2𝜂𝜁 ), uniformly in 𝑗, 𝑥 . This proves ( 30). At fixed 𝑥, the covariance identity for 𝑔 gives
E[ 𝐴2
𝑛


𝑉 𝑗 (𝑥)


2] ≤ 𝐶4− 𝑗. Summation over 𝑗 ≥ 𝐽𝑛 controls the supremum of the proxies. The common
error tends to zero on this tail, since 𝜀1/2−𝑑/4 ≤ 1 for 𝑑 ≤ 2. This proves ( 31), which is pointwise in 𝑥,
not a spatial supremum law. □
Lemma 4.3 (Residual comparison). Put 𝑄 𝑗 = 1
4

𝑔, 𝐶 𝐵2 𝑗 𝑔

, 𝑞 𝑗 = p𝑄 𝑗 and 𝐷 𝑛 = (𝑛𝜀3𝑑/2) −1/4. With
common 𝐸𝑛 = 𝑜P(1) and 𝑀𝑛 = 𝑂P(1), simultaneously for 𝑗 ≥ 2,
√𝑛
𝑟 𝑗 − 𝑞 𝑗
 ≤ 𝐸𝑛 + 𝑀𝑛{ 𝜀1/2−𝑑/42− 𝑗 + 𝐷 𝑛𝜀−𝑑/82− 𝑗 } . (33)
Proof. The exact mass field is 𝐹𝑗 = exp{( 𝑒 𝑗 + log 𝑃𝑛 (exp 𝑒 𝑗))/2} . Its linear part is (𝐼 + 𝑃𝑛)𝑒 𝑗/2, with
remainder bounded in sample 𝐿2 by 𝐶𝛿 𝑛


𝑒 𝑗



𝑛. Dividing by 𝑆 𝑗 centers this field, with the same remainder
order. Intertwining and (𝐼 + 𝐶)𝑅 = 𝐼 reduce the linear term to 𝑆∗
𝑛𝐵 𝑗 𝑔/2, with the common trajectory
error and 𝑂P(𝜂𝜁 ). Its mean obeys the exact identity

1, 𝑆∗
𝑛𝐵 𝑗 𝑔

𝑛 =

b𝑣 − 𝑣, 𝐵 𝑗 𝑔

= 𝑂P(𝜁 22− 𝑗),
because 𝑆𝑛1 = b𝑣 and 𝐵𝑣 = 0. It is root- 𝑛 negligible since √𝑛𝜁 2 = 𝜂 → 0. Finally


𝑆∗
𝑛ℎ



𝑛 − ∥𝑆∗ℎ∥ 𝐿2 ( 𝜇)
 ≤



b𝐶 − 𝐶




1/2
op
∥ℎ∥ H .
Apply this to 𝐵 𝑗 𝑔 and use √𝑛√𝜂 𝜁 = 𝐷 𝑛𝜀−𝑑/8. This proves ( 33). No uniformly root- 𝑛 negligible norm
comparison is asserted for the early iterates. □
8

5 Spectral concentration and the first hit
The proxy 𝑄 𝑗 is defined using the deterministic population operator, not a random spectral decomposition
of the sample matrix.
Lemma 5.1 (Hilbert quadratic concentration) . Let 𝜆𝑙 be the eigenvalues of 𝐶 and 𝛽𝑙 = (1 − 𝜆𝑙)/2. For
𝑗 ≥ 1, put 𝑇𝑗 = tr(𝐶2𝐵2 𝑗). Then
E𝑄 𝑗 = 𝑇𝑗
4𝑛 , Var𝑄 𝑗
(E𝑄 𝑗)2 ≤ 𝐶
max𝑙 𝜆2
𝑙 𝛽2 𝑗
𝑙
𝑇𝑗
+ 𝐶𝜀 −𝑑/24− 𝑗
𝑛𝑇𝑗
. (34)
The same calculation applies to deterministic spectral projections. On a band with 𝛽 ∈ [ 𝛽𝐻, 𝛽 𝐻, 2] ⊂
(0, 1/2) of fixed positive width, its relative variance is bounded by
𝐶 (𝜀𝑑/2 + 𝑛−1) (𝛽𝐻, 2/𝛽𝐻)2 𝑗 . (35)
Proof. For centered iid Hilbert vectors 𝑍𝑖 of covariance Σ, expansion of their squared sample-mean norm
gives
Var


𝑛−1 Í
𝑖 𝑍𝑖


2
= E ∥𝑍 ∥4 − ( tr Σ)2
𝑛3 + 2(𝑛 − 1) tr(Σ2)
𝑛3 . (36)
Only equal or reversed off-diagonal index pairs contribute. Terms with a singly occurring inde-
pendent centered variable vanish by conditioning; the diagonal term in ( 36) is retained. Take
𝑍 𝑗, 𝑦 = 1
2𝐶1/2𝐵 𝑗 [𝑎(𝑦)𝑘 𝑦 − 𝑣]. The stationary rank-one term vanishes for 𝑗 ≥ 1, so its covariance
is 𝐶2𝐵2 𝑗/4. The feature bound gives sup𝑦


𝑍 𝑗, 𝑦


2 ≤ 𝐶𝜀 −𝑑/24− 𝑗. Bound its fourth moment by this
envelope times its second moment, and use tr(Σ2) ≤ ∥Σ∥op tr Σ. This proves ( 34).
For a projected band, Lemma 3.2 gives 𝑇𝑗,𝐻 ≥ 𝑐𝜀 −𝑑/2 𝛽2 𝑗
𝐻 , with its 𝜆 values bounded away from zero.
The projected individual vector has the sharper envelope


𝑍 𝑗, 𝑦,𝐻


2 ≤ 𝐶𝜀 −𝑑/2 𝛽2 𝑗
𝐻, 2. In ( 36), this bounds
the diagonal relative term by 𝐶 (𝛽𝐻, 2/𝛽𝐻)2 𝑗/𝑛, and the covariance term by 𝐶𝜀 𝑑/2(𝛽𝐻, 2/𝛽𝐻)2 𝑗. This
proves (35). □
Proposition 5.2 (Stopping count). For fixed0 < 𝛾 < 𝑑 /4, the rule in (5) satisfies (7) and hits its threshold
strictly before 𝐿𝑛 with probability tending to one. In particular 𝜏𝛾 → ∞.
Proof. Put 𝑢∗ = 𝑑/4 − 𝛾 > 0. Fix a small slack 0 < 𝛿 < min(𝑢∗, 𝑑/8)/2, write 𝑢− = 𝑢∗ − 𝛿, and let
𝑗− =
 𝑢− log(1/𝜀)
log 2

, 𝑗 + =
 (𝑢∗ + 𝛿) log(1/𝜀)
log 2

.
Choose fixed 𝛽∗ < 1/2 close enough to 1/2 that, with 𝜎 = | log 𝛽∗| /log 2 > 1,
𝑢−𝜎 < 𝑢 ∗, 𝑢 − (𝜎 − 1) < 𝑑 /8. (37)
A fixed band with 𝛽 ≥ 𝛽∗ and 𝜆 bounded away from zero has 𝑐𝜀 −𝑑/2 eigenvalues, hence 𝑇𝑗 ≥ 𝑐𝜀 −𝑑/2 𝛽2 𝑗
∗ .
For every 2 ≤ 𝑗 ≤ 𝑗−, (34) is bounded by
𝐶 (𝜀𝑑/2 + 𝑛−1) (4𝛽2
∗) − 𝑗 .
The sum of these geometric bounds is 𝑜(1): its largest power is at most 𝐶𝜀 −2𝑢− ( 𝜎−1) , and ( 37) and
𝑛𝜀3𝑑/2 → ∞ control both terms. A Chebyshev union bound therefore gives, with probability tending to
one, simultaneously √𝑛 𝑞 𝑗 ≥ 𝑐𝜀 −𝑑/4 𝛽 𝑗
∗, 2 ≤ 𝑗 ≤ 𝑗−. (38)
The smallest lower bound, at 𝑗−, dominates 𝜀−𝛾, since 𝑢−𝜎 < 𝑢 ∗. The empirical-norm term in ( 33),
divided by (38), is bounded in probability by
𝐶 𝐷𝑛𝜀𝑑/8(2𝛽∗) − 𝑗− ≤ 𝐶 𝐷𝑛𝜀𝑑/8−𝑢− ( 𝜎−1) → 0.
9

The other root- 𝑛 errors are 𝑜P(1) or 𝑂P(1) for 𝑑 ≤ 2, whereas the smallest proxy lower bound diverges.
Thus no actual threshold hit occurs for 2 ≤ 𝑗 ≤ 𝑗−, and 𝑟2/𝑡𝑛 → ∞.
For the upper bracket use 𝑇𝑗 ≤ 4− 𝑗 tr(𝐶2) ≤ 𝐶𝜀 −𝑑/24− 𝑗. At 𝑗+, E[𝑛𝑄 𝑗+]/𝜀−2𝛾 ≤ 𝐶𝜀 2 𝛿 → 0, so
Markov’s inequality gives 𝑞 𝑗+/𝑡𝑛 → 0. The empirical-norm error divided by 𝑡𝑛 is 𝑂P(𝐷 𝑛𝜀𝑢∗+ 𝛿+𝛾−𝑑/8) =
𝑂P(𝐷 𝑛𝜀𝑑/8+ 𝛿) = 𝑜P(1). The other terms also vanish relative to 𝑡𝑛, proving a hit by 𝑗+. Moreover 𝑗+ < 𝐿 𝑛
eventually: the sampling condition in ( 2) bounds log(1/𝜀) above by 2 log 𝑛/(3𝑑) for large 𝑛.
It remains to exclude 𝑗 = 0, 1 without imposing a delay on the rule. For small log oscillation, replace
𝑒 by 𝑒 − 𝜇𝑛𝑒, which leaves the residual unchanged. Its linear centered mass field is (𝐼 + 𝑃𝑛)𝑒/2. On the
mean-zero subspace its norm lies between ∥𝑒∥𝑛 /2 and ∥𝑒∥𝑛. The nonlinear error is at most 𝐶 ∥𝑒∥∞ ∥𝑒∥𝑛.
Consequently there are fixed local constants such that 𝑟 (𝑒) is comparable to ∥𝑒∥𝑛 after centering, and
𝑟 (Ψ𝑛 (𝑒)) ≤ 𝐶𝑟 (𝑒). (39)
Indeed Ψ𝑛 kills constants, centering is an 𝐿2 contraction, and (24) controls its centered output. No nonlin-
ear mean-preservation or global residual monotonicity is assumed. Since ∥𝑒1∥∞ = 𝑜P(1) and 𝑟2/𝑡𝑛 → ∞,
(39) implies 𝑟1/𝑡𝑛 → ∞ and excludes 𝑗 = 1.
If the fixed 𝑓 is nonconstant, uniform consistency 𝐾𝑛1 → 𝑓 and the law of large numbers give
𝑟0 − →{ Var𝜇 (
p
𝑓 )} 1/2
𝜇
p
𝑓
> 0.
As 𝑡𝑛 → 0, this excludes 𝑗 = 0. If 𝑓 = 1, Lemma 3.3 gives ∥𝑒0∥∞ = 𝑜P(1), so (39) gives 𝑟0/𝑡𝑛 → ∞ and
excludes 𝑗 = 0 too. Thus the unrestricted first hit lies between 𝑗− and 𝑗+ with probability tending to one.
Letting the arbitrary slack 𝛿 tend to zero proves ( 7). □
Proposition 5.3 (Attained residual). For each fixed 0 < 𝑐 < 1/2, (8) holds.
Proof. For each realized 𝑔, 𝑄 𝑗 is a positive spectral sum. Its log-convexity implies that 𝑞 𝑗+1/𝑞 𝑗 is non-
decreasing and at most 1/2 when defined. Let 𝐽 = 𝑗− for the slack above. Fix 0 < 𝛽 𝐿 < 1/2, and choose
𝛽𝐿 < 𝛽 𝐻 < 𝛽 𝐻, 2 < 1/2 with sufficiently small positive width that
2𝑢− log(𝛽𝐻, 2/𝛽𝐻)/log 2 < 𝑑 /2.
The low- 𝛽 trace at 𝐽 is at most 𝐶𝜀 −𝑑/2 𝛽2𝐽
𝐿 , whereas the high band has trace at least 𝑐𝜀 −𝑑/2 𝛽2𝐽
𝐻 . The
latter concentrates by ( 35), including its diagonal term. The ratio of the two expectations tends to zero.
Markov’s inequality for the low part and concentration for the high part show that modes 𝛽 ≤ 𝛽𝐿 con-
tribute 𝑜P(𝑄 𝐽 ). Thus 𝑞 𝐽+1/𝑞 𝐽 ≥ 𝛽𝐿 − 𝑜P(1). Letting 𝛽𝐿 approach 1/2 and using monotonicity gives
sup
𝑗 ≥ 𝐽
𝑞 𝑗+1/𝑞 𝑗 − 1/2
 → 0 in probability. (40)
These ratios are defined on an event of probability tending to one, since the retained high band has positive
energy there.
Uniformly for 𝑗 ≥ 𝐽, ( 33) divided by 𝑡𝑛 is 𝑜P(1). Its empirical-norm term is bounded by
𝑂P(𝐷 𝑛𝜀𝑢−+𝛾−𝑑/8) = 𝑂P(𝐷 𝑛𝜀𝑑/8− 𝛿), and the commutator term is 𝑂P(𝜀1/2− 𝛿). At the first hit, 𝑟 𝜏𝛾 −1 > 𝑡 𝑛,
while (40) gives
𝑟 𝜏𝛾 ≥ ( 1/2 − 𝑜P(1))𝑡𝑛 − 𝑜P(𝑡𝑛).
The definition supplies 𝑟 𝜏𝛾 ≤ 𝑡𝑛, proving the claim. This does not assert a unique limit for 𝑟 𝜏𝛾 /𝑡𝑛; discrete
crossing may retain a rounding phase. □
6 Score limits and pointwise inference
Proposition 6.1 (Exact and stopped scores) . Under (2), at each fixed 𝑥 the exact score satisfies
𝐴𝑛{ b𝑠(𝑥) − ∇ log 𝑓 (𝑥)} ⇒ 𝑁𝑑 (0, [4𝐽𝑑/ 𝑓 (𝑥)] 𝐼𝑑).
The same limit holds at the stopped iterate in Theorem 2.1.
10

Proof. The reciprocal-branch expansion of Lemma 3.3 gives
b𝑠(𝑥) − 𝑠pop(𝑥) = 2 ⟨𝜃 𝑥, 𝑅𝑔⟩ + 𝑜P( 𝐴−1
𝑛 ), 𝑠 pop = −2∇ log 𝑎. (41)
Indeed b𝑠 = 2∇ logb𝑣 and the linear derivative is 2∇[ 𝑎(b𝑣 − 𝑣)] . The logarithmic nonlinear remainder is
bounded by 𝑂P(𝜂𝜀 −𝑑/4−1/2𝜁), and replacing b𝑣 − 𝑣 by 𝑅𝑔 has the same derivative bound. After multipli-
cation by 𝐴𝑛 both are 𝑂P(√𝑛𝜂𝜁 ) = 𝑜P(1).
The leading term in ( 41) is a centered iid triangular array. Its normalized covariance is
4𝜀𝑑/2+1 
𝜃 𝑥,𝜈 , 𝑅 (𝐶 − 𝑣 ⊗ 𝑣)𝑅𝜃 𝑥,ℓ

.
The rank-one subtraction is zero, because 𝑅𝑣 = 𝑣/2 and ⟨𝜃 𝑥, 𝑣⟩ = ∇[ 𝑎(𝑥)𝑣(𝑥)] = 0. Lemma 3.2 gives
the limiting covariance (4𝐽𝑑/ 𝑓 (𝑥)) 𝐼𝑑. A normalized summand is bounded in every fixed direction by
𝐶𝜀 𝑑/4+1/2
√𝑛 ∥𝜃 𝑥 ∥ H sup
𝑦


𝑎(𝑦)𝑘 𝑦 − 𝑣



H ≤ 𝐶√𝑛𝜀 𝑑/4 → 0.
The multivariate Lindeberg theorem therefore applies. Finally the population bias is 𝑂 (𝜀) and 𝐴𝑛𝜀 =√
𝑛𝜀 𝑑/2+3 → 0, proving the exact-score limit.
Proposition 5.2 places the stopping index beyond a deterministic sequence tending to infinity with
probability tending to one. Lemma 4.2 then gives ( 9) at that same-sample random index. Slutsky’s
theorem gives the stopped limit. □
Proof of Theorem 2.1. Propositions 5.2, 5.3 and 6.1 prove the count, attained-residual and score asser-
tions. The joint event in ( 8) follows by intersection of their events of probability tending to one. Its lower
bound also implies √𝑛 𝑟𝜏𝛾 → ∞. □
Corollary 6.2 (Feasible pointwise inference). Let e𝑓 (𝑥) = 𝐾𝑛1(𝑥) and let 𝑥 be fixed. Then
𝐴𝑛
q
e𝑓 (𝑥)
2√𝐽𝑑
{ 𝑠𝜏𝛾 (𝑥) − ∇ log 𝑓 (𝑥)} ⇒ 𝑁𝑑 (0, 𝐼𝑑), (42)
𝐴2
𝑛 e𝑓 (𝑥)
4𝐽𝑑


𝑠𝜏𝛾 (𝑥) − ∇ log 𝑓 (𝑥)


2 ⇒ 𝜒2
𝑑. (43)
Consequently the ball
8>>> <
>>>:
𝑢 ∈ R𝑑 :


𝑠𝜏𝛾 (𝑥) − 𝑢


 ≤
2
q
𝐽𝑑 𝜒2
𝑑, 1− 𝛼
𝐴𝑛
q
e𝑓 (𝑥)
9>>> =
>>>;
has limiting coverage 1 − 𝛼 for the fixed true score, for 0 < 𝛼 < 1. Predetermined directional intervals
follow from (42).
Proof. The periodic KDE is positive, its bias tends to zero, and its variance at fixed 𝑥 is 𝑂 (( 𝑛𝜀 𝑑/2) −1).
Hence e𝑓 (𝑥) → 𝑓 (𝑥) > 0 in probability. Slutsky’s theorem uses no independence between the plug-in
and the stopped score. The chi-square law and continuity at its quantile prove coverage. □
For 𝑑 = 2 the constant has the closed form𝐽2 = (𝜋2/12−log 2)/( 2𝜋): polar coordinates and 𝑢 = 2𝜋2𝑟 2
reduce it to (2𝜋) −1 ¯ ∞
0 𝑢(1 + 𝑒𝑢) −2 𝑑𝑢. Integration by parts, followed by the geometric-series integral for¯ ∞
0 𝑢/(1 + 𝑒𝑢) 𝑑𝑢, gives 𝜋2/12 − log 2.
11

7 Upper tolerance exponent
The coefficient of the threshold in ( 5) is one. The following statement concerns that convention and a
counterexample in the stated density class.
Theorem 7.1. Assume (2) and 𝑓 = 1. If 𝛾 ≥ 𝑑/4, then P(𝜏𝛾 = 0) → 1. At each fixed query,
𝐴𝑛𝑠𝜏𝛾 (𝑥) ⇒ 𝑁𝑑 (0, 𝑊𝑑), 𝑊 𝑑 = 2(4𝜋) −𝑑/2𝐼𝑑,
whereas 𝐴𝑛b𝑠(𝑥) ⇒ 𝑁𝑑 (0, 4𝐽𝑑 𝐼𝑑) and 𝑊𝑑 − 4𝐽𝑑 𝐼𝑑 is positive definite.
Proof. For uniform sampling, 𝑎 = 𝑣 = 1 and 𝐶 is the heat covariance operator. Write ℎ = 𝐾𝑛1 − 1 = 𝑔0.
The actual initial mass field is 𝐹0 =
√
1 + ℎ. Its linear part is ℎ/2. We have ∥ℎ∥∞ = 𝑂P(𝜂) and ∥ℎ∥ H =
𝑂P(𝜁), so its nonlinear 𝐿2 error is 𝑂P(𝜂𝜁 ). The empirical-to-population norm change is 𝑂P(√𝜂 𝜁 ), by
the sample heat covariance comparison. The empirical mean is exactly 𝜇𝑛ℎ = ∥ℎ∥2
H, because ⟨ℎ, 1⟩ = 0.
It is 𝑂P(𝜁 2) = 𝑜P(𝜁). Thus the actual centered and normalized residual satisfies
𝑟0 = 1
2 ∥ℎ∥ 𝐿2 (𝑑 𝑥) + 𝑜P(𝜁).
The squared leading norm has mean (4𝑛) −1 Í
𝑘≠0 𝜆2
𝑘. Its relative variance is 𝑂 (𝜀𝑑/2): in ( 36) the indi-
vidual squared feature norm is constant by translation invariance, so its diagonal variance vanishes, andÍ 𝜆4
𝑘/(Í
𝑘≠0 𝜆2
𝑘)2 = 𝑂 (𝜀𝑑/2). Gaussian lattice sums consequently yield
𝑟0
𝑛−1/2𝜀−𝑑/4 − →𝑐𝑑 = 1
2 (4𝜋) −𝑑/4 < 1 in probability. (44)
For 𝛾 ≥ 𝑑/4 the threshold is at least 𝑛−1/2𝜀−𝑑/4, so the actual rule stops at zero with probability tending
to one.
Here 𝑠0 = 2∇ log(1 + ℎ). At a fixed query, 𝐴𝑛∇ℎ is tight by its kernel second moment, and ℎ(𝑥) → 0.
Its kernel triangular-array CLT has covariance
𝑊𝑑 = 4
¹
R𝑑
(2𝜋𝜉 ) (2𝜋𝜉 )𝑇 𝑒−4 𝜋2 | 𝜉 | 2
𝑑𝜉 = 2(4𝜋) −𝑑/2𝐼𝑑.
The Lindeberg bound is again 𝐶/(√𝑛𝜀 𝑑/4) → 0; the logarithmic denominator is negligible by Slutsky.
The exact-score limit follows from Proposition 6.1. Writing Λ(𝜉) = 𝑒−2 𝜋2 | 𝜉 | 2
, its covariance is
4𝐽𝑑 𝐼𝑑 = 4
¹
(2𝜋𝜉 ) (2𝜋𝜉 )𝑇 [Λ(𝜉)/( 1 + Λ(𝜉))] 2 𝑑𝜉.
For every nonzero direction, the difference from 𝑊𝑑 is the integral of a squared nonzero linear projection
against a strictly positive radial weight. It is positive, proving the claim. □
The theorem does not classify arbitrary endpoint multipliers or assert failure for every nonconstant
density. It establishes the upper exponent obstruction for a uniform guarantee over the model class with
this rule.
8 Scope of the inference statements
The confidence regions above are pointwise-in- 𝑓 asymptotic regions at predetermined spatial queries.
They do not provide uniform finite-sample coverage over the unrestricted positive 𝐶6 class. A standard
discontinuous-functional obstruction [9] makes that distinction necessary. For example, at a query 𝑥0 and
each fixed real 𝑞,
𝑓𝑚,𝑞 (𝑥) = 1 + 𝑞
2𝜋𝑚 sin(2𝜋𝑚 (𝑥1 − 𝑥0,1))
12

has true score 𝑞𝑒 1 at 𝑥0, is smooth and lies between 1/2 and 3/2 for large integer 𝑚, while its 𝑛-sample
total variation distance from uniform data is at most 𝑛| 𝑞|/( 2𝜋2𝑚). For each fixed 𝑛 this tends to zero.
Any uniformly honest confidence interval would therefore have to cover every fixed𝑞 under uniform data
with its nominal probability; an almost surely bounded interval cannot do so. Quantitative smoothness in-
formation or a specified smoothed target is needed for such finite-sample claims. This does not contradict
Theorem 2.1, which fixes the density before taking its limit.
The pointwise result permits confidence regions for a density’s logarithmic gradient and tests of a
zero gradient at a prespecified point. It does not locate an unknown mode or justify data-selected queries.
Uniform spatial limits, boundaries, vanishing densities and diffusion-path error control require additional
arguments. At 𝑑 = 3 the two conditions on 𝑛𝜀9/2 in ( 2) are incompatible; higher dimensions are not
covered. The leading update law is an asymptotic iteration count, and finite-sample runtime improvements
are not inferred from it.
AI assistance
AI tools assisted preparation and checking. The author is responsible for the content.
References
[1] J. Feydy, T. Séjourné, F.-X. Vialard, S.-i. Amari, A. Trouvé and G. Peyré. Interpolating between
optimal transport and MMD using Sinkhorn divergences. In Proceedings of the 22nd International
Conference on Artificial Intelligence and Statistics, Proceedings of Machine Learning Research 89
(2019), 2681–2690. https://proceedings.mlr.press/v89/feydy19a.html.
[2] A. Kokot. Variational Problems in Feature Selection and Data Compression . Ph.D. dissertation,
University of Washington, 2026. https://hdl.handle.net/1773/57893.
[3] G. Mordant. The statistics of entropic optimal transport with decreasing regularisation. Preprint,
arXiv:2412.12007v2, 29 April 2026. https://arxiv.org/pdf/2412.12007v2.
[4] C. L. Wormell and S. Reich. Spectral convergence of diffusion maps: improved error bounds and
an alternative normalization. SIAM Journal on Numerical Analysis 59(3) (2021), 1687–1734.
doi:10.1137/20M1344093.
[5] X. Cheng and B. Landa. Bi-stochastically normalized graph Laplacian: convergence to manifold
Laplacian and robustness to outlier noise. Information and Inference: A Journal of the IMA 13(4)
(2024), iaae026. doi:10.1093/imaiai/iaae026.
[6] A. Celisse and M. Wahl. Analyzing the discrepancy principle for kernelized spectral filter learning
algorithms. Journal of Machine Learning Research 22(76) (2021), 1–59.
https://jmlr.org/papers/v22/20-358.html.
[7] I. Maeda, R. Y ao and A. Nitanda. Statistical analysis of the Sinkhorn iterations for two-sample
Schrödinger bridge estimation. Preprint, arXiv:2510.22560v1, 26 October 2025.
https://arxiv.org/abs/2510.22560v1.
[8] S. Schubr. Residual certificates and empirical sharpness for entropic self-transport scores.
Luminary Research Archive, edition 1, 18 October 2026. Permanent archive.
[9] J.-M. Dufour. Some impossibility theorems in econometrics with applications to structural and
dynamic models. Econometrica 65(6) (1997), 1365–1388. doi:10.2307/2171740.
13
