NONUNIFORM SELF-SINKHORN STOPPING IN A TRUE-SCORE INFERENCE WINDOW
Research proof candidate, 18 October 2026. Independent review required.
All earlier proofs, code, experiments and publication packages are read only.

1. Exact theorem

Let f be a fixed positive C^6 probability density on the unit circle T,
with min f>0, and let X_1,...,X_n be iid with law mu=f dx. The density may
be nonconstant. Use the exact normalized periodic Gaussian k_epsilon with
Fourier multipliers exp(-2 pi^2 epsilon k^2). Assume
  epsilon=epsilon_n ->0,
  n epsilon^(3/2) ->infinity,
  n epsilon^(7/2) ->0.                                      (H)
For epsilon=n^-kappa this gives the nonempty interval 2/7<kappa<2/3.
No optimal bandwidth assertion is made.

Let a_hat be the exact empirical self-scale, extended off sample by
a_hat K_n a_hat=1. Run the actual damped self-Sinkhorn iteration
  b_0=1, b_(j+1)=sqrt(b_j/(K_n b_j)),
where K_n integrates against mu_n=n^-1 sum delta_Xi. Let
  s_j=2 d_x log(K_n b_j), s_hat=2 d_x log(K_n a_hat),
  F_j=sqrt(b_j K_n b_j), S_j=int F_j dmu_n,
  r_j=||F_j/S_j-1||_(L^2(mu_n)).
Scores and residuals are invariant under positive rescaling of b_j.

Two updates are mandatory. Then stop at the first j>=2 with
  r_j<=t_n, t_n=n^-1/2 epsilon^-1/16.
To define the procedure on every sample, impose the budget
L_n=max(3,ceil(4 log n)); if no hit by L_n, return iteration L_n.
Call this index tau_n. The two-update requirement is part of this theorem;
no first-hit assertion from j=0 or j=1 is being substituted for it.

Set A_n=sqrt(n) epsilon^(3/4) and
  J_n=floor[log(1/epsilon)/(16 log 2)].
For each fixed x in T, the following hold:
(a) With probability tending to one, J_n<tau_n<L_n, the threshold is
    actually hit, and t_n/8<=r_(tau_n)<=t_n. In particular
       sqrt(n) r_(tau_n) ->infinity in probability.
(b) A_n |s_(tau_n)(x)-s_hat(x)| ->0 in probability.
(c) Both exact and stopped scores estimate the TRUE nonconstant score:
       A_n [s_hat(x)-f'(x)/f(x)] => N(0,4 J_1/f(x)),
       A_n [s_(tau_n)(x)-f'(x)/f(x)] => N(0,4 J_1/f(x)),       (1)
    where J_1=int_R (2 pi xi)^2/[1+exp(2 pi^2 xi^2)]^2 dxi>0.

This is a one-dimensional periodic theorem in a genuine undersmoothing
window. It does not invoke a Euclidean boundary-model CLT. The increment
over the earlier uniform theorem is fixed nonuniform density, control of
population burn-in, and a window in which population score bias is negligible.

2. Population facts and a common feature space

Write H for heat convolution on Lebesgue measure, a for the population
self-scale, v=1/a=K_mu a, and w=Hf=K_mu 1. Standard heat-kernel estimates
and a local contraction about a=f^-1/2 give
  ||a sqrt(f)-1||_(C^2)<=C epsilon.                          (2)
In particular a,v and their derivatives through order two are uniformly
bounded, with a,v bounded away from zero. The same assertion in C^1 gives
  ||-2(log a)'-(log f)'||_infinity<=C epsilon.                 (3)

For completeness, (2) uses q=a sqrt(f), g=sqrt(f), and the equation
q H(qg)=g in C^2. The multiplier of (I+H)^-1 is 1-F(sqrt(epsilon)k),
where F(xi)=1/[1+exp(2 pi^2 xi^2)] is Schwartz. Periodizing its inverse
Fourier transform bounds this inverse uniformly on C^2. The derivative
of q H(qg) at q=1 is g(I+H)+E_epsilon, with ||E_epsilon||_(C^2->C^2)
=O(sqrt(epsilon)); commute each derivative through H and use the Gaussian
first absolute moment on each commutator. The forcing H g-g is O(epsilon)
in C^2 for g in C^4. The remaining term is quadratic, so a fixed small
ball contraction gives (2). C^6 is more than sufficient. Uniqueness follows
from PSD: for two positive solutions a,b,
  <a-b,K_mu(a-b)>=-int (a-b)^2/(ab) dmu<=0.

Let H_epsilon be the real RKHS of k_epsilon. A feature k_y has norm
sqrt(K_0)=O(epsilon^-1/4); a derivative feature has norm
sqrt(K_1)=O(epsilon^-3/4). Let
  S u=int a(y) k_y u(y) dmu(y), S* h(y)=a(y)h(y),
  C=S S*=int a(y)^2 k_y tensor k_y dmu(y), P=S* S.
P is PSD, self-adjoint and Markov on L^2(mu); hence ||S||=||C||=1.
Also ||v||_H=1, Cv=v and S*v=1. Put
  B=(I-C)/2, R=(I+C)^-1.
Then ||B||<=1/2 and ||R||<=1. Define the derivative section
  theta_x=d_x[a(x)k_x], ||theta_x||_H<=C epsilon^-3/4.

We shall use the population derivative Gram limit
  epsilon^(3/2) <theta_x,R C R theta_x>_H ->J_1/f(x).           (4)
Here is its reduction. Under Uu=sqrt(f)u, P is unitarily conjugate to
Q=M_q H M_q with q=a sqrt(f)=1+O(epsilon) in C^1. Thus ||Q-H||=O(epsilon)
and their I-plus resolvents differ by O(epsilon). The transformed derivative
row differs from f(x)^-1/2 d_x k(x,.) in L^2(dx) by
O(sqrt(epsilon) epsilon^-3/4). Fourier orthogonality then reduces (4) to
epsilon^3/2 sum_k (2 pi k)^2[lambda_k/(1+lambda_k)]^2 ->J_1.
This is a Gaussian-dominated lattice Riemann sum. This is also the population
Gram lemma in the preserved corrected manuscript; no empirical conclusion
from that manuscript is needed as a new premise here.

The nonzero eigenvalues of C and P agree, and hence those of C and Q agree.
Consequently the ordered eigenvalues of C differ from those of H by
O(epsilon), by the compact self-adjoint min-max principle. Every fixed
positive-width spectral interval strictly inside (0,1), counted using a
smaller interior band of heat eigenvalues, therefore contains at least
c epsilon^-1/2 eigenvalues for small epsilon. Moreover
  tr C=K_0 int a^2 dmu<=C epsilon^-1/2.                        (5)
Only these fixed-band counts, not a changing-frequency Weyl expansion, are
used below.

3. Quantitative empirical scaling in the RKHS

All random-order bounds in the rest of the proof concern this same sample.
They will hold on a common event with probability tending to one, uniformly
in iteration j. Put the deterministic rates
  zeta=n^-1/2 epsilon^-1/4, eta=n^-1/2 epsilon^-1/2,
and define centered feature averages
  g_n=(mu_n-mu)(a k), g_0=(mu_n-mu)(k).
Second moments give ||g_n||_H+||g_0||_H=O_p(zeta). In fact
  Cov(g_n)=(C-v tensor v)/n, <v,g_n>_H=0 exactly.               (6)
Let C_0n=int a^2 k tensor k dmu_n. The Hilbert--Schmidt second moment
bound gives
  ||C_0n-C||_op<=||C_0n-C||_HS=O_p(eta),                     (7)
because E||a(X)^2 k_X tensor k_X||_HS^2<=C K_0^2. No logarithmic
concentration bound or independent-query assumption is needed for (7).

Consider F_n(z)=z-K_n(1/z) near v as a map on the RKHS. At v,
F_n(v)=-g_n and D F_n(v)=I+C_0n, whose inverse has norm at most one.
For ||h||_H<=C||g_n||_H, the embedding bound gives ||h||_infinity=O_p(eta).
As eta->0, v+h stays positive. Expanding the reciprocal gives a nonlinear
remainder of H-norm at most C sqrt(K_0)||h||_H^2, with the corresponding
Lipschitz bound on that ball. To see this bound, write the remainder as a
feature integral with factor h^2/[v^2(v+h)]; one h is bounded in sup norm,
and the other is controlled by the empirical covariance C_0n. Its operator
norm is O_p(1) by (7). Thus a contraction with fixed linear inverse proves
the local positive solution and
  h_n:=v_hat-v=R g_n+O_p(eta zeta) in H,
  ||h_n||_H=O_p(zeta),                                      (8)
where v_hat=K_n a_hat. The branch is the exact empirical scale: its site
values solve the unique positive finite scaling equations, and the equation
fixes the same off-sample extension. No global signed-measure claim is used.
The replacement of (I+C_0n)^-1 by R costs O_p(eta zeta) by (7).

Embedding (8) and differentiating the feature section give
  ||a_hat/a-1||_infinity=O_p(eta),
  ||a_hat-a||_infinity+sqrt(epsilon)||(a_hat-a)'||_infinity
       =O_p(eta).                                           (9)
Define exact empirical maps
  S_n u=int a_hat(y)k_y u(y)dmu_n(y),
  S_n* h(y)=a_hat(y)h(y) at sample sites,
  C_hat=S_n S_n*, P_n=S_n* S_n.
P_n is PSD and Markov, so ||S_n||=||C_hat||=1. Equations (7)--(9),
or the quadratic-form comparison of the two sets of feature weights, give
  ||C_hat-C||_op=O_p(eta).                                   (10)
The section theta_hat_x=d_x[a_hat(x)k_x] satisfies, uniformly in x,
  ||theta_hat_x||_H<=O_p(1)epsilon^-3/4,
  ||theta_hat_x-theta_x||_H=O_p(eta epsilon^-3/4).             (11)
All norms of evaluation of H-functions in L^2(mu_n) are uniformly bounded
by their H-norm on this event, using min a>0 and ||C_0n||=O_p(1).

4. The first update, its commutator, and population burn-in

Let e_j=log(b_j/a_hat), including the off-sample extension. The exact update
is
  e_(j+1)=[e_j-log(P_n exp(e_j))]/2.                          (12)
The first update is b_1=(K_n 1)^-1/2. Put
  d_1=log v-(1/2)log w.
From (2), ||d_1||_(C^2)=O(epsilon). Using (8), the logarithm expansion,
and L^2(mu_n) evaluation bounds,
  e_1=d_1+a R g_n-(1/2)g_0/w+O_(L^2(mu_n),p)(eta zeta).

The commutator q_n=g_0-v g_n has kernel expression
  q_n(x)=(mu_n-mu)_y k(x,y)[1-v(x)a(y)].
The bracket vanishes at y=x and is O(dist(x,y)), uniformly. Given X_i=x,
the sample's diagonal term vanishes EXACTLY. Its conditional mean is
  -[w(x)-v(x)^2]/n=O(epsilon/n),
and its conditional variance is at most
  (C/n) int k(x,y)^2 dist(x,y)^2 dmu(y)<=C epsilon^1/2/n.
Therefore
  ||q_n||_(L^2(mu_n))=O_p(sqrt(epsilon) zeta).
This accounts explicitly for dependence between query sites and the sample.
Also v/w-a=O(epsilon). Since R-(1/2)I=B R, (9) now gives
  e_1=d_1+S_n* B R g_n+u_1,
  ||u_1||_(L^2(mu_n))=O_p[(eta+sqrt(epsilon))zeta],
  ||e_1||_(L^2(mu_n))=O_p(epsilon+zeta),
  ||e_1||_infinity=O_p(epsilon+eta).                          (13)

Let B_n=(I-P_n)/2. The population operator satisfies
  ||(I-P)d_1||_infinity<=C epsilon ||d_1||_(C^2)=O(epsilon^2).
This follows by Taylor expansion of the smooth test function in a normalized
Gaussian row, including its smooth first-moment drift. In the empirical row,
changing a to a_hat costs O_p(eta)||d_1||_infinity=O_p(epsilon eta).
The feature average for a d_1 has H-norm fluctuation O_p(epsilon zeta),
and for a has O_p(zeta). Embedding bounds and the row denominator bounded
away from zero therefore give
  ||B_n d_1||_(L^2(mu_n))=O_p(epsilon^2+epsilon eta).          (14)
This second update removes the remaining deterministic population transient
at the statistical accuracy used below.

5. An actual all-iteration expansion after two updates

The Markov property in (12) gives the exact sup-norm nonexpansion
||e_(j+1)||_infinity<=osc(e_j)/2<=||e_j||_infinity. Write
delta_n=||e_1||_infinity=O_p(epsilon+eta)=o_p(1).
For ||e||_infinity<=delta_n, Taylor expansion in the empirical L^2 space,
Jensen's inequality and Markov contraction show
  ||Psi_n(e)-B_n e||_2<=C delta_n ||e||_2.
As P_n is PSD, ||B_n||_2<=1/2. On one event with probability tending to
one the entire trajectory is therefore L^2-contractive with a fixed factor
q=3/4. Summing its nonlinear recursion proves, uniformly j>=1,
  ||e_j-B_n^(j-1)e_1||_2
       <=O_p[(epsilon+eta)(epsilon+zeta)].                    (15)
This uses separate L^2 and sup control, rather than incorrectly inferring
relative-direction control from a potential norm alone.

The feature intertwining identity is exact:
  B_n S_n*=S_n* B_hat, B_hat=(I-C_hat)/2.
Both ||B_hat|| and ||B|| are at most 1/2. Telescoping their powers and (10)
give a uniform error O_p(eta zeta) when applied to B R g_n. Combining
(13)--(15), and using (14) for j>=2, yields
  e_j=S_n* B^j R g_n+E_(j,n), j>=2,
  sup_(j>=2)||E_(j,n)||_(L^2(mu_n))=O_p(rho_n),              (16)
where one may take
  rho_n=epsilon^2+epsilon eta+epsilon zeta+eta zeta
                              +sqrt(epsilon) zeta.
Assumption (H) gives sqrt(n)rho_n->0: its five terms after multiplication
by sqrt(n) are respectively
  sqrt(n)epsilon^2, epsilon^1/2, epsilon^3/4,
  (n epsilon^3/2)^-1/2, epsilon^1/4.
The first tends to zero because n epsilon^7/2->0. All bounds in (16) use
one common event and apply simultaneously to every iteration j>=2.

6. Point-score and residual consequences of the expansion

The exact score error is 2 d_x log(P_n exp(e_j)). The derivative row has
L^2(mu_n) norm at most ||theta_hat_x||_H. Expanding the logarithm at small
e_j, its difference from 2 d_x(P_n e_j) is bounded by
C epsilon^-3/4 delta_n ||e_j||_2. Equations (10)--(11) and (16) therefore
give, uniformly in j>=2 and spatial x,
  s_j(x)-s_hat(x)=V_j(x)+o_p(A_n^-1),
  V_j(x)=2<theta_x,C B^j R g_n>_H.                           (17)
The operator/section replacements cost O_p(epsilon^-3/4 eta zeta);
after multiplication by A_n their size is (n epsilon^3/2)^-1/2=o(1).

For the residual, expand the square-root mass field in L^2(mu_n):
  F_j=1+(I+P_n)e_j/2+O_2(delta_n ||e_j||_2).
Using (16) and (10), its linear term is
  (1/2)S_n* B^j g_n+o_(L^2(mu_n),p)(n^-1/2),
uniformly j>=2. The mean of this linear term is particularly simple:
  <1,S_n* B^j g_n>_n=<v_hat,B^j g_n>_H
       =<v_hat-v,B^j g_n>_H=O_p(zeta^2 2^-j),
because Bv=0 and j>=1. It is o_p(n^-1/2). Thus centering and dividing
by S_j preserve the stated residual-norm error. Write
  Q_j=(1/4)<g_n,C B^(2j)g_n>_H.
Comparing the empirical and population quadratic forms by (10) gives the
KEY bound, simultaneously j>=2:
  |r_j-sqrt(Q_j)| <= o_p(n^-1/2)
                       +O_p(sqrt(eta) zeta 2^-j).             (18)
The random constants are common to all j. This is deliberately a norm bound.
It does NOT assert a root-n o_p(1) residual approximation uniformly in j;
that stronger claim would be false as a deduction from the estimates here.

For the score proxy, (6) gives at a fixed query x
  E[A_n^2 |V_j(x)|^2]<=C4^-j,
since ||C||,||R||<=1 and ||theta_x||<=Cepsilon^-3/4. Therefore
  E sup_(j>=J) A_n^2 |V_j(x)|^2<=C4^-J.
Equation (17) implies
  sup_(j>=max(2,J_n)) A_n |s_j(x)-s_hat(x)| ->0
in probability. Spatial uniformity is claimed only for the remainder in
(17); the last stochastic limit is at a FIXED query.

7. Spectral residual concentration with nonuniform observations

Let lambda_l be the eigenvalues of C, beta_l=(1-lambda_l)/2, and j>=1.
Equation (6) and Bv=0 imply
  E Q_j=T_j/(4n), T_j=tr(C^2 B^(2j)).                        (19)
The Perron eigenvalue contributes zero for j>=1.
This trace formula does not assert independent empirical coordinates.

Here is the variance bound needed instead. Define iid centered H-vectors
  Z_(j,y)=(1/2)C^1/2 B^j[a(y)k_y-v].
Then Q_j=||n^-1 sum_i Z_(j,Xi)||_H^2. For any iid centered Hilbert vector Z
with covariance Sigma,
  Var||n^-1 sum Z_i||^2
   =[E||Z||^4-(tr Sigma)^2]/n^3
          +2(n-1)tr(Sigma^2)/n^3.                           (20)
This follows by separating diagonal terms and equal/reversed off-diagonal
index pairs. It does not omit the variable diagonal present for nonuniform f.
Here Sigma_j=(1/4)C^2 B^(2j), and
  sup_y ||Z_(j,y)||^2<=C epsilon^-1/2 4^-j.
Consequently
  Var Q_j/(E Q_j)^2
     <=C max_l(lambda_l^2 beta_l^(2j))/T_j
                  +C epsilon^-1/2 4^-j/(n T_j).              (21)
The same proof applies after any deterministic spectral projection of C.

By (5), the fixed band lambda_l in [1/4,1/2] has at least c epsilon^-1/2
eigenvalues, and beta_l>=1/4 there. Thus, with J=J_n,
  T_J>=c epsilon^-1/2 16^-J>=c epsilon^-1/4,
  max_l lambda_l^2 beta_l^(2J)<=4^-J<=Cepsilon^1/8.
The relative variance in (21) is O(epsilon^3/8+1/(n epsilon^1/8))=o(1).
It follows that
  4n Q_J/T_J ->1 in probability,
  sqrt(n Q_J)>=c epsilon^-1/8 with probability tending to one. (22)

For each realized g_n, Q_j is a positive spectral sum. It is decreasing
and log-convex; hence Q_(j+1)/Q_j is nondecreasing when Q_j>0.
To bound that ratio below after J, split at beta_l=1/4. The low trace
is at most Cepsilon^-1/2 16^-J. A high band lambda_l in [1/8,1/4] has
trace at least c epsilon^-1/2(9/64)^J, and its relative variance tends to
zero by (21) restricted to that band: the bounds include
  Cepsilon^1/2(7/6)^(2J) +C(7/6)^(2J)/n ->0.
The low/high expected ratio is at most C(4/9)^J->0. Markov for the low
part and concentration for the high band show that the low contribution
is o_p(Q_J). Thus, with probability tending to one,
  Q_(j+1)>=Q_j/32 for EVERY j>=J.                            (23)
This argument requires no diagonalization of an empirical random operator.

8. Actual first-hit control

Set D_n=(n epsilon^3/2)^-1/4->0. The extra measure term in (18) has
root-n size O_p(n^-1/4 epsilon^-1/2 2^-j).
For 2<=j<=J, bound it by its value at j=2 or simply drop 2^-j. Relative
to the lower bound epsilon^-1/8 in (22), its size is O_p(D_n)=o_p(1).
The o_p(1) root-n remainder is also negligible relative to that bound.
Since Q_j>=Q_J, (18) excludes every threshold crossing in [2,J]
with probability tending to one; epsilon^-1/8 dominates sqrt(n)t_n
=epsilon^-1/16.

For j>=J, retain 2^-j<=Cepsilon^1/16. Dividing (18) by t_n gives
  sup_(j>=J)|r_j-sqrt(Q_j)|/t_n ->0 in probability,            (24)
because n^-1/4 epsilon^-1/2 epsilon^1/16 epsilon^1/16=D_n.
This is a threshold-relative approximation, which is exactly what the
stopping proof needs throughout the full bandwidth window (H).

The L^2 contraction following (15) bounds r_j by
O_p(epsilon+zeta)(3/4)^(j-1). At L_n-1 this is o_p(n^-1/2), so the
threshold is hit strictly before the budget with probability tending to
one. Also J_n<L_n eventually by n epsilon^3/2->infinity.
At the first genuine hit tau>J, r_(tau-1)>t_n. Equations (23)--(24) give
  r_tau >=(1/sqrt(32))(1-o_p(1))t_n-o_p(t_n)>=t_n/8
with probability tending to one, while r_tau<=t_n by definition. This
proves (a). The fixed-query simultaneous score tail from (17) evaluated
at this same random tau proves (b); no optional-stopping independence is used.

9. True-score CLT in the same nonuniform periodic model

From (8), the exact score fluctuation is
  s_hat(x)-s_pop(x)=2<theta_x,R g_n>_H+o_p(A_n^-1),
where s_pop=-2(log a)'. Its leading term is a centered iid triangular
array with covariance
  4/n <theta_x,R(C-v tensor v)R theta_x>
       =4/n <theta_x,R C R theta_x>.
Indeed Rv=v/2 and <theta_x,v>=d_x[a(x)v(x)]=0. By (4), the scaled
variance tends to 4J_1/f(x). The influence function is bounded by
C||theta_x|| sqrt(K_0)<=Cepsilon^-1. Thus each normalized summand is
bounded by C/[sqrt(n)epsilon^1/4]->0, verifying Lindeberg. The triangular-
array CLT yields the centered Gaussian limit for the exact empirical score.
By (3), its population bias is O(epsilon), and
  A_n epsilon=sqrt(n epsilon^7/2)->0.
This proves the first limit in (1); (b) and Slutsky prove the second.
The true score is not assumed zero. This completes the theorem.

10. What is and is not established

Prior public work already contains the two-update population mechanism:
Alex Kokot, Variational Problems in Feature Selection and Data Compression,
University of Washington dissertation, publicly deposited 16 September 2026,
Section 4.4.1, Lemma 4.4 and Corollary 4.3. Stable primary record:
https://hdl.handle.net/1773/57893 . Under q=b^-2, its symmetric update is
the same population half-self update, up to the Gaussian parameter convention.
Chapter 3, Lemma 3.29 also explicitly uses a Hilbert CLT and the diagonal
(I+T)^-1 scaling derivative at fixed regularization. Neither that method nor
the population two-update cancellation is claimed as a new invention here.

The proof supplies fixed nonuniform density, an explicit useful inference
window, actual standard-start iterates with two mandatory updates, and the
attained residual lower bound at an observable first hit. In contrast to
the preceding uniform proof, (18) is NOT uniformly root-n accurate: its
decaying measurement term must be kept through the stopping argument.

Known tools include scaling, RKHS covariance calculus, implicit equations,
the half-step spectral filter, kernel CLTs, spectral min-max, Hilbert-vector
quadratic variance, log-convexity and Slutsky. Originality of the theorem's
combination and the value of its precise conclusion require separate review.
No claim is made for dimensions >=2, a Euclidean support with boundary,
optimal tolerance, optimal bandwidth, measured runtime improvement, or
external human peer review. The 18 finite experiments are diagnostic only;
this analytic theorem does not use their floating outputs as premises.
