Sinkhorn Linearization and the Spectral Proxy: Unifying the Statistical and Algorithmic Theory of Feature-Parameterized Inverse Optimal Transport via a Single Spectral Sandwich

Sinkhorn Linearization and the Spectral Proxy: Unifying the Statistical and Algorithmic Theory of Feature-Parameterized Inverse Optimal Transport via a Single Spectral Sandwich Han Dong1,† , Jiaming Li2,† , Yongqiang Gong1 , Ruixi Li1 , Yin Liu1 1 School of Medicine, Nankai University 2 College of Artificial Intelligence, Nankai University

†These authors contributed equally to this work. dh411424@163.com, xmkevin2004@gmail.com, gongyq@mail.nankai.edu.cn, 2211996@mail.nankai.edu.cn, liuyin tiger@outlook.com

arXiv:2608.13201v1 [stat.ML] 13 Aug 2026

August 14, 2026 Abstract We develop the statistical and algorithmic theory of inverse optimal transport (IOT) under the feature-parameterized cost Cθ (i, j) = −θ ⊤ φ(i, j), where the observations are the conditional transition operators of the entropic optimal transport plan. The core technical contribution is the Sinkhorn linearization—the implicit-function sensitivity of the entropic OT plan to the cost, obtained by differentiating the KKT conditions—together with its spectral proxy, a formula that is spectrally exact (it preserves all singular-value bounds) yet geometrically transparent. The restricted Hessian of entropic OT on the tangent space satisfies the spectral sandwich −1 (πmin /ε)I ⪯ HTp ⪯ (πmax /ε)I, from which the single core spectral bound σmin (Jθ ) ≥ (πmin /(amax ε)) λmin (Σ) is derived, driving the entire theory. On this core we establish four theorems and one observation. T1 (identifiability): θ is identifiable on the quotient of the gauge kernel (RF /NΦ ), with the dimension bound F ≤ (K − 1)2 and the rank condition rank(Σ) = F ; global injectivity is proved rigorously by a three-step composition argument. T2 (sparsistency): the ℓ1 -penalized estimator recovers the true support under all-coordinate score concentration, irrepresentability of the actual Hessian, and a global selection condition; the Hoeffding exponent for a single normalized empirical average is 2nt2n /∆2max , and with multiple marginals the rate is re-scaled by a weighted effective scale. T3 (well-posedness): we introduce the feature-moment map M (θ) = Φ⊤ xθ and proceed in two parts: (a) local strong monotonicity is derived from the core spectral bound without compactness; (b) global strong monotonicity additionally requires a compact parameter domain; the Lipschitz bound of the Q-space inverse map is L ≤ ε∥Φ⊤ Sa ∥op /(πmin λmin (Σ)). T4 (convergence): local strong convexity of the cross-entropy objective is derived from the core spectral bound, with pa2 rameter µ ≥ πmin λmin (Σ)/ε2 ; gradient descent converges monotonically to a local minimum; the empirical Hessian concentrates at the population Hessian at rate O(n−1/2 ) under regularity bounds. O5 (misspecification): under a compact parameter domain, a uniform law of large numbers, and a unique pseudo-true point, the estimator converges to the projection of the truth onto the OT model set; the Hölder continuity of the projection map remains a conjecture, and numerics report only a setting-dependent empirical effective exponent αeff . Dependency structure. T1 (identifiability) is the foundation; T3 and T4 are both corollaries of the core spectral bound (Section 3), forming a well-posedness-and-convergence package; T2 requires the additional IR condition (T1 + IR); O5 uses T1 + T3. All identifiability and consistency statements are conditional on a fixed ε; identifiability is taken modulo the scale p coupling (θ, ε) 7→ (cθ, cε). The effective information scale also depends on πmin (θ, ε) λmin (Σ)/ε; ε ≫ cond(Σ)−1 is only an empirical heuristic, not a threshold rigorously derived in this paper.

Key words. inverse optimal transport; Sinkhorn linearization; spectral proxy; identifiability; sparsistency; well-posedness; convergence; misspecification; entropic regularization 1

AMS subject classifications. 49Q22, 62F12, 62J07, 90C25

1

Introduction

The problem. Optimal transport (OT) computes a mass-transfer plan given a cost and two marginals; inverse optimal transport (IOT) goes the other way, recovering the cost from observed transport plans. In applications such as biological lineage tracing, transfer learning, and economic matching, what is observed is typically the induced conditional transition operator rather than an explicit cost, and one must invert the data back to the cost parameters that drive the transitions. This paper treats the feature-parameterized cost case, in which the cost is a linear combination of interpretable features. This formulation turns IOT into a statistical inference problem: recover a sparse cost parameter from finitely many samples of conditional transition operators. Core contribution: Sinkhorn linearization and the spectral proxy. The central technical contribution of this paper is the Sinkhorn linearization (Section 3)—the implicitfunction sensitivity of the entropic OT plan to the cost. Implicit differentiation of the KKT conditions yields the exact linear response of the plan to a cost perturbation, δx = −BHT−1 B ⊤ δc, where the restricted Hessian HT satisfies the spectral sandwich (πmin /ε)I ⪯ HT−1 ⪯ (πmax /ε)I. From this we introduce the spectral proxy δxSSP = −(1/ε)PT Dπ PT δc—a formula that keeps the same upper and lower spectral bounds given by πmin , πmax while remaining geometrically transparent, reading as “project to the tangent space → multiply elementwise by π → p project back to the tangent space”. The core spectral bound σmin (Jθ ) ≥ (πmin /(amax ε)) λmin (Σ) follows directly from the spectral sandwich and drives all of the theorems T1–T4 and Observation O5. Overview of the method. We formulate IOT as a cross-entropy maximum-likelihood problem on entropic transport plans and establish four theorems and one observation. Identifiability (T1) delineates in what sense the parameter can be uniquely recovered; sparsistency (T2) gives sufficient conditions for ℓ1 -penalized estimation together with an explicit exponential rate constant; well-posedness (T3) introduces the feature-moment map M (θ) = Φ⊤ xθ and establishes the stability of the inverse map in two parts; convergence (T4) analyzes the local convergence dynamics of the optimization; and the misspecification analysis (O5) answers where the estimator goes when the data-generating mechanism departs from the OT assumption. Difference from existing work. Existing IOT work is mostly formulated with an explicit or low-dimensional cost, and typically lacks explicit characterizations of the geometric degeneracy that identifiability owes to the gauge subspace, of the exponential rate constant for sparse recovery, and of the convergence target of the estimator under misspecification. The contribution of this paper is to reduce these ingredients to computable constants: the featuredimension bound F ≤ (K − 1)2 for identifiability (a prior-structural constant), the exponential rate constant C2 = 2/∆2max (with t2n kept in the exponent; a posterior-diagnostic constant), the Lipschitz constant LΘ ≤ ε∥Φ⊤ Sa ∥op /(πmin λmin (Σ)) (derived from the Sinkhorn linearization, with compactness guaranteeing the global bound; a posterior-diagnostic constant) together with its numerical estimate, and the comparison of the misspecification projection residual against random models. Comparison with classical IOT work. [24] first posed the IOT problem and gave a Bayesian inference framework; the present paper builds a frequentist statistical theory on top of that foundation. [25] propose Bayesian IOT: a prior is placed on the cost and posterior consistency is derived under a correctly specified model; their identifiability relies on a Gaussian-process prior on the cost, whereas our T1 gives a finite-dimensional rank condition; their consistency is Bayesian (posterior concentration), ours is frequentist (ℓ1 support recovery with an exponential rate, T2); we give an explicit rate constant, which they do not address. [13] study statistical inference for nonparametric OT maps estimated by kernel density esti2

mation; their setting is forward (estimate the map given the cost), ours is inverse (recover the cost given the map), and their sample-complexity analysis does not address the gaugedegeneracy obstruction at the core of T1. [12] develop center-outward distribution functions via OT, providing a notion of multivariate quantiles; their identifiability is trivial (the map is defined by construction), whereas our T1 must contend with the gauge non-uniqueness that arises when recovering the cost from the plan. The economics line ([8], [9]) establishes IOT identifiability under parametric utility models but addresses neither sparse recovery (T2) nor misspecification (O5); our O5 applies directly to economic settings in which observed matching data are almost certainly not OT-generated. Recent IOT progress. [11, 10] systematically study identifiability conditions for nonlinear and discrete IOT, giving necessary-and-sufficient linear-programming combinatorial characterizations; our T1 supplements these with an explicit dimension bound and rank condition in the feature-parameterized setting. [1] establish irrepresentability conditions for sparse recovery in IOT and bridge to the graphical Lasso; our T2 gives an explicit exponential rate constant within a unified spectral framework. [2] analyze the well-posedness and algorithms of Bregman-regularized IOT; our T3 focuses on the spectralproxy-driven derivation of the Lipschitz bound. [18] apply IOT to semi-supervised learning; [15] propose a convex-optimization form of inverse estimation for Markov chains. These works emphasize different aspects; the distinctive contribution of the present paper is the unified spectral framework provided by the Sinkhorn linearization, from which all four theorems and one observation are derived from a single core spectral bound. Strength of the claims. The four theorems and one observation have different degrees of completeness in their arguments. T1 (the gauge structure and the global injectivity argument), the local part of T3 (local strong monotonicity, derived from the core spectral bound), and T4 (Hessian positive definiteness derived from the core spectral bound, with 7.2 bridging to finite samples) are conclusions with complete derivation chains. T2 relies on irrepresentability of the actual OT information matrix, all-coordinate score concentration, and a global selection condition; Lemma 5.2 gives an a priori verifiable condition on the feature Gram matrix, but its transfer to the OT Hessian requires an additional bridge. The global part of T3 holds on a compact parameter domain (Assumption 6.6 is required). The entropic bias of T3 satisfies the local order O(|ε′ − ε|), derived from differentiability (6.3) and the implicit function theorem, under a fixed parameter scale (excluding the joint scale coupling) and an invertibility assumption on the cross-entropy Hessian, with zero bias at ε′ = ε. The convergence of O5 is given by the misspecified M -estimation theory of [28]; the Hölder continuity of the projection map is a conjecture. The identifiability of ε is given by 2.1, and πmin has an ε-dependent lower-bound characterization. The open problems are collected at the end of this introduction. Open problems. The following are the known limitations. Conjectures: (i) Hölder continuity of the projection map Q 7→ arg minQ′ ∈MOT CE(Q∥Q′ ); the αeff in the experiments is only a finite-setting empirical effective exponent, not a universal theoretical exponent; (ii) random OT models are typically far from the projection, so that the residual bound of O5 is non-vacuous. Open: (iii) global strong monotonicity on the unbounded parameter space RF /NΦ (without compactness); (iv) the expression of the optimal exponential rate constant in terms of cond(Σ), the sparsity s, and the irrepresentability margin; (v) the geometry of the basins of attraction of global optimization in T4. The experiments are in Section 9. Notation. θ ∈ RF is the cost parameter; φ(i, j) ∈ RF is the feature vector; K is the number of states; ε is the entropic regularization level; G is the gauge subspace of row and column constants; x = vec(π) is the vectorized plan; c = vec(C) is the vectorized cost.

2

Preliminaries

Entropic optimal transport. Given a cost matrix C ∈ RK×K , marginals a, b ∈ ∆K (all entries strictly positive: ai > 0, bj > 0), and a regularization parameter ε > 0, the entropic OT 3

plan is π(C, a, b) = arg

min

n o X ⟨C, π⟩ + ε πij log πij ,

π∈M(a,b)

(1)

ij

where M(a, b) = {π ≥ 0 : π1 = a, π ⊤ 1 = b} is the set of couplings with marginals a, b, with the convention 0 log 0 := 0. The unique solution has the Sinkhorn form π = diag(u) exp(−C/ε)diag(v) ([23]; [7]). Note that the scaling vectors u, v carry a scalar gauge redundancy: (u, v) 7→ (cu, c−1 v) leaves the plan unchanged. We work with the feature-parameterized cost Cθ (i, j) = −θ ⊤ φ(i, j). For comprehensive treatments of OT see [19], [22]; for the connection between entropic OT and Schrödinger problems see [14]. Observation model. The observations are the conditional transition operators Qi→j = πij /ai . For each marginal pair (al , bl ) there are n independent observations, giving the empirical (n) b (n) = Nl,ij /Nl,i (with the convention joint plan Pbl and the empirical conditional operator Q l,ij b (n) is filled with the prior uniform distribution that when Nl,i = 0 the corresponding row of Q l 1/K; for n sufficiently large the probability of Nl,i = 0 decays exponentially and does not affect the asymptotic conclusions). The marginals b al , bbl of the empirical joint plan generally differ from the prescribed marginals al , bl ; in the theoretical analysis of this paper the prescribed marginals are treated as known. The conditional transition operator is the actually observable quantity, also referred to as state-transition data. The IOT estimation problem. Given L marginal pairs (al , bl ) and observed conditional b l , the IOT estimator minimizes the cross-entropy with an ℓ1 penalty: operators Q min

θ∈RF

L X

b l ∥ Qθ (al , bl ) + λn ∥θ∥1 , CE Q

(2)

l=1

where Qθ (al , bl ) := Q(Cθ , al , bl ) is the conditional transition operator induced by the parameterized cost, and λn is a penalty parameter scaled with the P sample P bsize. The cross-entropy b between conditional operators is defined as CE(Q∥Q) = − i ρi j Q ij log Qij ; in this paper we take ρi = 1 (equal row weights), consistent with the population objective of Theorem 7.1. If the samples are drawn at the joint frequencies, the natural choice ρi = ai (source-marginal weights) is also possible. The gauge subspace and the identifiability obstruction. If a function of the form ui +vj (a sum of row and column constants) is added to the cost, the Sinkhorn plan is unchanged. Hence the cost is identifiable only modulo G; at the parameter level, the object of identification is the gauge kernel NΦ , obtained by pulling G back to the parameter space through Φ (see Section 4). G is defined as n o dim G = 2K − 1. (3) G := C ∈ RK×K : C = 1α⊤ + β1⊤ , 2

Let G ∈ RK ×(2K−1) be an orthonormal basis of G and ΠG = GG⊤ the orthogonal projection onto G (acting on vectorized matrices, with the Frobenius inner product ⟨A, B⟩F = tr(A⊤ B)). ⊤ The gauge-cleaned feature vectors are φ⊥ k := (I − GG ) vec(Φk ), and the gauge-cleaned Gram matrix is h iF ⊥ Σ := φ⊥ , φ ∈ RF ×F . (4) ′ k k ′ k,k =1

⊥ If NΦ ̸= {0}, we write λ+ min (Σ) for the smallest positive eigenvalue of Σ on NΦ ; under Assumption 4.3, λ+ min (Σ) = λmin (Σ).

Proposition 2.1 (Identifiability of ε under normalization). For (θ, ε) ∈ RF ×R>0 , identifiability holds up to the joint rescaling (θ, ε) 7→ (cθ, cε). The equivalence relation is Qθ,ε = Qθ′ ,ε′ ⇐⇒ 4

θ θ′ − ′ ∈ NΦ , ε ε

where NΦ is the gauge kernel (Section 4). Therefore, under a normalization constraint (e.g. ∥θ∥ = 1), ε is uniquely identifiable from the conditional transition operator Qθ,ε , provided NΦ = {0} (i.e. the columns of the feature matrix Φ are linearly independent after modding out G). If NΦ ̸= {0}, unit norm does not remove the non-uniqueness along the gauge-kernel directions: there exist θ ̸= θ ′ with ∥θ∥ = ∥θ ′ ∥ = 1 and θ/ε−θ ′ /ε′ ∈ NΦ , for which Qθ,ε = Qθ′ ,ε′ . Hence the identifiability of ε must be distinguished: it is identifiable on the quotient space RF /NΦ , and identifiable on the original space RF only when NΦ = {0}. There is an additional global scale coupling beyond G: the plan depends only on C/ε = −θ ⊤ φ/ε, so (θ, ε) 7→ (cθ, cε) leaves the plan unchanged; identifiability statements are made conditional on a fixed ε, and the identifiability of ε itself and of its joint scale with θ is resolved by Proposition 2.1 (ε is uniquely identifiable under a normalization constraint). Remark 2.2 (A uniform lower bound for the Sinkhorn plan entries). The constant πmin := minij πij appears in the Lipschitz bound L ≤ ε/(πmin λ+ min (Σ)) (the T3 bound on the quotient + space; in the full-rank case λmin (Σ) = λmin (Σ)) and in the score-concentration bound of T2. For a single marginal pair and a single feature coordinate, (3.6) gives ∆k ≤ 2

πmax ∥PT φk ∥2 . επmin

In the multi-marginal case, ∆max should be defined according to the sampling weights, the marginal indices, and the coordinate set in T2. We now give the correct limiting behavior of πmin as ε → 0 together with an exponential-type lower bound. Basic setup. Let a, b ∈ ∆◦K (all entries strictly positive), C ∈ RK×K , ε > 0. Write h(π) = P i,j πij log πij ; then πε = arg min {⟨C, π⟩ + εh(π)}. π∈M(a,b)

Since the Gibbs kernel e−C/ε is strictly positive entrywise, πε,ij > 0 and the optimal solution is unique. Entropy-selected limit. Let the optimal value and the optimal set of the unregularized OT be C ⋆ = min ⟨C, π⟩, U = arg min ⟨C, π⟩. π∈M(a,b)

π∈M(a,b)

Define the entropy-selected solution on U: π H = arg min h(π), π∈U

unique by the strict convexity of h. For any π̄ ∈ U, the optimality of πε gives ⟨C, πε ⟩ + εh(πε ) ≤ C ⋆ + εh(π̄). Since ⟨C, πε ⟩ ≥ C ⋆ ,   0 ≤ ⟨C, πε ⟩ − C ⋆ ≤ ε h(π̄) − h(πε ) . h is bounded on the transport polytope (− log(K 2 ) ≤ h(π) ≤ 0), so 0 ≤ ⟨C, πε ⟩ − C ⋆ ≤ ε log(K 2 ), and every limit point belongs to U. Meanwhile h(πε ) ≤ h(π̄); passing to the limit shows that all limit points minimize h on U. Since this minimizer is unique, πε −→ π H

(ε ↓ 0).

5

If the unregularized OT solution is unique, π H is that unique solution; otherwise the limit is the max-entropy solution on the optimal face. Exponential lower bound. Let ∆C = maxi,j Cij − mini,j Cij , Rε = e∆C /ε . The P Sinkhorn form is πε,ij = ui Kij vj with Kij = e−Cij /ε . Since Rε−1 ≤ Kij /Kij ′ ≤ Rε , writing Sj = i ui Kij , the column marginal bj = vj Sj gives vj = bj /Sj , hence   bj Sj ′ bj vj −1 bj = ∈ Rε , Rε . vj ′ bj ′ Sj bj ′ bj ′ Therefore

  Kij vj πε,ij −2 bj 2 bj = ∈ Rε , Rε . πε,ij ′ Kij ′ vj ′ bj ′ bj ′

In particular, πε,ij ′ ≤ Rε2 (bj ′ /bj )πε,ij . Summing over j ′ and using Rε2 πε,ij /bj , i.e. πε,ij ≥ ai bj e−2∆C /ε ,

P

j ′ πε,ij ′

= ai gives ai ≤

πmin (ε) ≥ amin bmin e−2∆C /ε .

Why a polynomial lower bound cannot hold. Take  
0 ∆ a = b = 21 , 12 , C = . ∆ 0
1 1 −∆/ε 1 e−∆/ε , so π A direct solution gives πε = 2(1+e1−∆/ε ) e−∆/ε . For any min (ε) = 2(1+e∆/ε ) ∼ 2 e 1 α c > 0 and finite α ≥ 0, πmin (ε)/(cε ) → 0, so no uniform bound of the form πmin (ε)
≥ cεα exists. This remains so even when the cost entries are all distinct (e.g. C = 02 14 , with πmin (ε) ≍ e−1/(2ε) ). Parameterized cost. As θ varies,   2 osc(Cθ ) πmin (θ, ε) ≥ amin bmin exp − , ε

osc(Cθ ) = max Cθ,ij − min Cθ,ij . i,j

i,j

For uniform use over a parameter domain Θ, one needs supθ∈Θ osc(Cθ ) < ∞. Summary. For fixed strictly positive marginals and a finite cost matrix, the entropic plan πε is strictly positive for ε > 0 and converges as ε ↓ 0 to the unique max-entropy solution π H on the optimal face. πmin satisfies an exponential-type lower bound and cannot be replaced by a uniform polynomial bound cεα . If the limiting plan contains zero entries, then πmin (ε) → 0, at a rate that depends on the optimal face, the dual gap, and the combinatorial structure of the cost matrix. Parameter domain and the local/global hierarchy. The cost parameter θ lives on the identifiable quotient space RF /NΦ (NΦ is the gauge kernel; see Section 4), which is unbounded. The core spectral bound (Section 3, (3.3)) is derived pointwise (it holds at each θ individually). On any compact convex subset the pointwise bound can be uniformized and yields strong monotonicity; on the unbounded quotient space the pointwise infimum may vanish and a global uniform bound may fail. T3 is stated in two parts: (a) local strong monotonicity (derived by uniformizing the pointwise bound on any compact convex subset); (b) global strong monotonicity (valid after restricting θ to a compact convex subset Θ ⊂ P RF /NΦ ). b l ∥Qθ ) weights Row weighting of the cross-entropy objective. The objective l CE(Q the transition distribution of each source state i equally (the operator is treated as K row distributions), which differs from the log-likelihood of n i.i.d. transition samples (weighted by the source frequencies ai ); this paper adopts the operator-level objective in order to treat all states symmetrically. 6

3

Sinkhorn Linearization and the Spectral Proxy

This section is the technical core of the paper. It derives the Sinkhorn linearization—the implicit-function sensitivity of the entropic OT plan to the cost—and introduces the spectral proxy, a formula that preserves the same upper and lower spectral bounds given by πmin , πmax while remaining geometrically transparent. All subsequent theorems (T1–T4, O5) build on the quantities defined here. For the analysis of the penalization trajectory of entropic OT see [27] and [6]; for the stability of entropic OT plans see [4]. [3] independently study the second-order curvature of continuous OT with respect to the cost and apply it to IOT identifiability; the discrete Sinkhorn linearization here is complementary to their work.

3.1

KKT optimality and implicit differentiation

Fix marginals a, b > 0 and ε > 0. The entropic OT objective is L(π) = ⟨C, π⟩ + ε

X

πij (log πij − 1),

(5)

ij

whose gradient is ∇π L = C + ε log π (entrywise logarithm). The KKT conditions for the constraint π ∈ M(a, b) are Cij + ε log πij + αi + βj = 0, where the αi + βj span the gauge subspace G. Equivalently, projecting onto the tangent space T = G ⊥ :
PT C + ε log π = 0 .

(6)

Equation (6) defines π as an implicit function of C. Differentiating along a perturbation δC (in vectorized form x = vec(π), c = vec(C)):
PT δc + εDπ−1 δx = 0,

δx ∈ T ,

(7)

where Dπ = diag(vec(π)). Interior smoothness (the implicit function theorem). To remove the one-dimensional redundancy among the row and column marginal constraints, the implicit function theorem 2 argument below takes A ∈ R(2K−1)×K to be the marginal-constraint matrix with one redundant constraint deleted, which has full row rank. The KKT conditions of the constrained problem are c + ε log x + A⊤ λ = 0, Ax = r, and its derivative matrix with respect to (x, λ) is  −1  εDπ A⊤ K(x) = . (8) A 0 If K(x)(u, v)⊤ = 0, then Au = 0 and u⊤ (εDπ−1 u + A⊤ v) = εu⊤ Dπ−1 u = 0 (because A⊤ v is orthogonal to ker A). Since Dπ−1 ≻ 0 we get u = 0, and then v = 0 because A has full row rank. Hence K(x) is nonsingular, and the implicit function theorem guarantees that, at fixed positive marginals and fixed ε > 0, x is a smooth function of the cost c.

3.2

The restricted Hessian and the spectral sandwich

The unconstrained Hessian ∇2x L = εDπ−1 is diagonal (the entropy is elementwise separable). But restriction to the tangent space changes the inverse structure: in general (P AP )|−1 T ̸= P A−1 P for a projection P . 2 2 Let B ∈ RK ×(K−1) be an orthonormal basis of T (B ⊤ B = I, range(B) = T ). Every feasible perturbation is δx = Bδz. The restricted Hessian in tangent-space coordinates is 7

HT = εB ⊤ Dπ−1 B ≻ 0.

(9)

Although HT is generally dense (not diagonal), its inverse satisfies a spectral sandwich, which is the key to all subsequent analysis. We state it as a lemma. 2

2

Lemma 3.1 (Spectral sandwich for the restricted Hessian). Let B ∈ RK ×(K−1) be an orthonormal basis of the tangent space T , and HT = εB ⊤ Dπ−1 B the restricted Hessian in tangentspace coordinates. Then its inverse satisfies πmin πmax I ⪯ HT−1 ⪯ I, ε ε

πmin := min πij , πmax := max πij . ij

ij

(10)

2

Proof. For any unit vector z ∈ R(K−1) , the Rayleigh quotient of HT = εB ⊤ Dπ−1 B at z is z ⊤ HT z (Bz)⊤ Dπ−1 (Bz) =ε , ⊤ z z (Bz)⊤ (Bz) where we used z ⊤ z = (Bz)⊤ (Bz), which holds because B ⊤ B = I and hence ∥Bz∥ = ∥z∥. Since Bz ∈ T and every entry of the probability matrix π satisfies 1/πij ∈ [1/πmax , 1/πmin ], the weighted average (Bz)⊤ Dπ−1 (Bz)/∥Bz∥2 lies in [1/πmax , 1/πmin ]. Therefore every eigenvalue λ of HT satisfies ε πmax

∥z∥2 ≤ z ⊤ HT z ≤

ε πmin

∥z∥2 ,

and inverting the matrix reverses the order of its eigenvalues (each eigenvalue of HT−1 is the reciprocal of one of HT ), so πmin /ε ≤ λ(HT−1 ) ≤ πmax /ε for every eigenvalue. This is exactly the matrix inequality (10). Note that HT is generally dense (not diagonal), and its inverse HT−1 = ε−1 (B ⊤ Dπ−1 B)−1 ̸= This is precisely where the tangent-space restriction changes the “project-theninvert-elementwise” structure. ε−1 B ⊤ Dπ B.

3.3

Exact linearization and the spectral proxy

From (7) in tangent-space coordinates: δx = −BHT−1 B ⊤ δc .

(11)

This is the exact Sinkhorn linearization. Its equivalent Schur-complement form (obtained by eliminating the Lagrange multipliers from the full KKT system) is δx = −

i 1h Dπ − Dπ A⊤ (ADπ A⊤ )−1 ADπ δc, ε

(12)

where A is the reduced constraint matrix of full row rank defined above, so that ADπ A⊤ is invertible. If instead one uses the matrix Afull containing all 2K row and column marginal constraints, its rank is 2K − 1, and the corresponding Schur complement should replace the + ordinary inverse by the Moore–Penrose pseudoinverse (Afull Dπ A⊤ full ) . Throughout this paper, all formulas use the reduced A and the ordinary inverse. The condition number of the reduced matrix may degenerate as πmin → 0; whether the limit is singular depends on whether the constraint graph of the positive-support entries remains sufficiently connected. Schur-complement derivation. Differentiating the KKT system c + ε log x + A⊤ λ = 0 along δc gives εDπ−1 δx + A⊤ δλ = −δc and Aδx = 0. Multiplying the first equation by ε−1 Dπ and left-multiplying by A, the constraint Aδx = 0 yields ADπ A⊤ δλ = −ADπ δc; solving for δλ and substituting A⊤ δλ back gives the Schur-complement form (12). 8

The spectral proxy (SSP). Although (11) is exact, it lacks the geometric transparency of a simple formula. We therefore introduce the Sinkhorn spectral proxy: 1 δxSSP := − PT Dπ PT δc . ε

(13)

Definition 3.2 (Spectral proxy). Formula (13) is a spectral proxy for the exact linearization (11): it is not elementwise equal to the exact derivative, but it is spectrally equivalent—for any δc, πmin πmax ∥PT δc∥ ≤ ∥δx∥ ≤ ∥PT δc∥, ε ε

πmin πmax ∥PT δc∥ ≤ ∥δxSSP ∥ ≤ ∥PT δc∥. ε ε

(11) and (13) share exactly the same spectral sandwich (10). Consequently, every statistical conclusion that depends only on spectral quantities (singular-value bounds, Lipschitz constants, strong-convexity parameters, rate constants) is identical under the exact linearization and under the spectral proxy. The spectral proxy has three desirable properties that the exact formula lacks: 1. Geometric interpretability: it reads as “project to the tangent space → multiply elementwise by π → project back to the tangent space”—each step has an independent geometric meaning. 2. Computational simplicity: evaluating the proxy requires only two tangent-space projections and an elementwise multiplication, with no (K − 1)2 × (K − 1)2 linear-system solve. 3. Intuition preservation: it reveals that the local sensitivity of the Sinkhorn map is driven by the plan entries themselves—large entries amplify the response, small entries suppress it. The exact formula (11) (or its Schur-complement form (12)) is used when elementwise derivatives are required (e.g. for computing the actual score in gradient descent). The spectral proxy is used for all operator-norm analyses and for conveying geometric intuition.

3.4

Feature parameterization and the Jacobian

P Under the cost parameterization Cθ = − Fk=1 θk Φk , we have δc = −Φ δθ, where Φ = [vec(Φ1 ), . . . , vec(ΦF )] ∈ 2 RK ×F . The conditional transition operator is Qij = πij /ai , vectorized as q = Ra x with −1 Ra = diag(a−1 1 , . . . , aK ) ⊗ IK . Throughout this paper we adopt the row-major (row-stacking) vectorization convention: vec(M ) = (M11 , M12 , . . . , M1K , M21 , . . . , MKK )⊤ . Under this convention Ra = diag(a−1 ) ⊗ IK correctly maps vec(π) to vec(Q); with column-major vectorization the Kronecker order would reverse to IK ⊗ diag(a−1 ). The Jacobian Jθ := ∂q/∂θ can be written in exact or spectral-proxy form: 1 JθSSP = Ra PT Dπ PT Φ. ε

Jθexact = Ra BHT−1 B ⊤ Φ,

Derivation. Define the plan sensitivity operator Sπ := BHT−1 B ⊤ . It is symmetric, positive semidefinite, and range(Sπ ) = vec(T ). From (11), ∂x/∂θ = Sπ Φ, and composing with q = Ra x gives Jθ = Ra Sπ Φ. Since the smallest singular value of Ra is σmin (Ra ) = 1/amax (amax := maxi ai ), the spectral sandwich (10) propagates to a singular-value lower bound on Jθ (see the singular-value bound below).

9

3.5

Singular-value lower bound (the core spectral result)

From the spectral sandwich (10) and the gauge-cleaned Gram Σ = Φ⊤ PT Φ, we obtain the core spectral bound that drives the entire theory: Proposition 3.3 (Singular-value lower bound). At every θ with πij (θ) > 0 and λmin (Σ) > 0, σmin (Jθ ) ≥

πmin (θ) p λmin (Σ) , amax ε

σmax (Jθ ) ≤

πmax (θ) p λmax (Σ) . amin ε

Proof. For any ∆θ ∈ RF , let δc = −Φ∆θ. From the spectral sandwich (10) and the exact formula (11): ∥δx∥ = ∥BHT−1 B ⊤ δc∥ ≥

πmin ⊤ πmin πmin p λmin (Σ) ∥∆θ∥. ∥B δc∥ = ∥PT Φ∆θ∥ ≥ ε ε ε

Since q = Ra x and σmin (Ra ) = 1/amax , we have ∥Jθ ∆θ∥ ≥ (1/amax )∥δx∥, giving the lower bound. The upper bound is analogous. Since amax ≤ 1, a simpler but looser bound is p σmin (Jθ ) ≥ (πmin /ε) λmin (Σ).

3.6

Score bound

The per-sample score at observation (i, j) is sk (i, j) = −∂θk log πij = −(∂θk πij )/πij . From the exact formula and the spectral sandwich (10): |sk (i, j)| ≤

πmax πmax ∥PT φk ∥2 ≤ ∥φk ∥2 , ε πmin ε πmin

where φk = vec(Φk ). This bound is the absolute-value envelope Mk of the score. The per-sample score has mean zero under correct specification, so its range (needed by Hoeffding) is ∆k := max sk − min sk ; because the score takes values on both sides (sk ∈ [−Mk , Mk ]) we have ∆k ≤ 2Mk , and the range must not be underestimated as half the envelope. This bound is more conservative than the naive elementwise formula (it carries an extra 1/ε factor and the ratio πmax /πmin ), but it is rigorously derivable from the spectral analysis of the restricted Hessian. When πmin and ε are bounded away from zero in a regular region, this gives a uniform score bound ∆max , which feeds into the Hoeffding/Bernstein concentration of T2. A sharper vector expression is available from the Schur-complement form (12): i 1h I − A⊤ (ADπ A⊤ )−1 ADπ φk , ε where gk = ∂x/∂θk . This reveals that the score is a Dπ -weighted gauge-corrected feature residual, not a simple Frobenius projection residual. Dπ−1 gk =

3.7

From pointwise to uniform: the conditioning hierarchy

The bound (3.3) is pointwise: it holds at each θ individually. The statistical theory distinguishes four levels of conditioning: 1. Pointwise identifiability (T1). At each θ with Σ ≻ 0 and πij > 0, the Jacobian has full column rank. No compactness is needed. 2. Local conditioning (T3a, T4). On any compact neighborhood K, inf θ∈K πmin (θ) > 0 is guaranteed by continuity, giving the uniform bound σmin ≥ τK /(amax ε) on K. This is a derived property, not an assumption.

10

3. Global uniform conditioning (T3b). This requires the entire parameter domain Θ to be compact with inf Θ πmin > 0 (Assumption 6.6). Without compactness the infimum may vanish even though every pointwise bound is positive. 4. Boundary degeneration (ε → 0 or πmin → 0). As the plan concentrates on the sparse OT support, πmin → 0and the  spectral bound degenerates. A 2 × 2 example: take 0 ∆ a = b = (1/2, 1/2), C = (∆ > 0); then ∆ 0 πmin (ε) =

1 1 ∼ e−∆/ε , ∆/ε 2 2(1 + e )

ε ↓ 0,

so πmin (ε)/ε → 0 exponentially. This shows that pointwise identifiability at a fixed regularization level is not equivalent to uniform statistical estimability as ε ↓ 0; this is the regime in which IOT statistical inference becomes ill-conditioned.

3.8

Relation to the downstream theory

The following table summarizes how each theorem depends on the quantities defined in this section: Theorem

Depends on

From this section

T1 (identifiability) T2 (sparsistency) T3 (well-posedness) T4 (convergence) O5 (misspecification)

ker(J ) = NΦ , rank(Σ) = F Score bound ∆max , rate constant C2 Lower bound on σmin (J ) λmin (∇2 ℓ) via σmin (J ) T1 + T3 for the regular region R

(6), (3.3) (3.6) (3.3) (3.3) (3.3), (10)

Table 1: Dependence of each theorem on the definitions of this section All five theoretical components radiate from the single spectral bound (3.3), which itself follows from the spectral sandwich (10) of the restricted Hessian. The spectral proxy (13) provides an equivalent, geometrically transparent formula that gives the same spectral bounds without solving a (K − 1)2 × (K − 1)2 linear system.

4

T1: Identifiability

4.1

Statement

Theorem 4.1 (Identifiability). Fix a regularization level ε > 0. Under the feature-parameterized cost Cθ , assume the marginal pairs are strictly positive and nondegenerate (Assumption 4.2: al,i > 0, bl,j > 0 for all i, j, and both sides have full support). Identifiability is stated in two layers: (0) Quotient-space identifiability (no full rank needed). The map [θ] ∈ RF /NΦ 7→ Qθ is always globally injective on the quotient space: Qθ = Qθ′ ⇐⇒ θ − θ ′ ∈ NΦ . This injectivity does not require rank(Σ) = F . (1) Gauge invariance (no full rank needed). If the cost difference falls into the gauge subspace, i.e. Φ(∆θ) ∈ vec(G), then the forward map is unchanged: Qθ+∆θ = Qθ . (2) Original-space identifiability (if and only if full rank). On the original parameter space RF , θ 7→ Qθ is injective if and only if rank(Σ) = F (equivalently NΦ = {0}). This full-rank condition only says that the parameter itself is uniquely identifiable in the original parameter space; it does not alter the injectivity of (0) on the quotient space.

11

Proof of global injectivity: by a three-step composition argument—(1) the linear parameterization θ 7→ Cθ is injective modulo the gauge; (2) the Sinkhorn map C 7→ π(C, a, b) is injective modulo the gauge (strict convexity of the dual → unique potentials → unique plan); (3) the normalization π 7→ Q is linear injective (ai > 0, guaranteed by Assumption 4.2). The composition of three injective maps is globally injective. Identification structure and the gauge kernel. The genuinely unidentifiable directions in the parameter are characterized by the parameter-induced gauge kernel  NΦ := v ∈ RF : Φv ∈ vec(G) , i.e. those parameter directions along which the cost falls into the gauge subspace G. The natural object of identification for the parameter is the quotient space RF /NΦ ; the object of identification for the cost is RK×K /G. The gauge-cleaned Gram is Σ = Φ⊤ PT Φ = (B ⊤ Φ)⊤ (B ⊤ Φ), and v ⊤ Σv = ∥PT Φv∥2 ,

Σ ≻ 0 ⇐⇒ NΦ = {0} ⇐⇒ rank(Σ) = F.

Since rank(PT ) = (K − 1)2 , Σ ≻ 0 implies the feature dimension bound F ≤ (K − 1)2 . Cost-level injectivity (Step 2 expanded). If C, C ′ produce the same plan π, then by the KKT conditions there exist multipliers (α, β) and (α′ , β ′ ) with Cij + ε log πij + αi + βj = 0 ′ + ε log π + α′ + β ′ = 0. Subtracting the two equations shows that C − C ′ has the and Cij ij i j form −αi P − βj , hence C − C ′ ∈ G. Conversely, for Gij = ui + vj and any π ∈ M(a, b), P ⟨G, π⟩ = i ui ai + j vj bj is independent of π, so adding G to the objective only adds a constant and does not change the unique optimal plan. Thus C 7→ π(C; a, b) is globally injective on RK×K /G. Parameter-level identification (Steps 1 and 3 merged). Since ai > 0, Qθ = Qθ′ ⇐⇒ πθ = πθ′ . By cost-level injectivity, Qθ = Qθ′ ⇐⇒ Cθ − Cθ′ ∈ G ⇐⇒ Φ(θ − θ′ ) ∈ vec(G) ⇐⇒ θ − θ′ ∈ NΦ . Hence θ is globally identifiable on RF /NΦ ; when NΦ = {0} (i.e. Σ ≻ 0), θ is uniquely identifiable on the original space RF . The effect of theP number of marginals on conditioning. The kernel of the stacked (l)⊤ J (l) is identically N (independent of L), so adding marginals sensitivity Gram SL = L Φ l=1 J does not change the identification rank. Let V := NΦ⊥ , and write λ+ min (Σ) for the smallest positive eigenvalue of Σ on V. For each marginal define the spectral envelopes on the quotient space     πmin,l 2 + πmax,l 2 + + ml := Ml := λmin (Σ), λmax (Σ). al,max ε al,min ε ! ! PL L L X X Ml+ + + ml IV ⪯ SL V ⪯ M l IV , condV (SL ) ≤ Pl=1 . L + l=1 ml l=1 l=1 The above condition number is defined on V (or equivalently on the quotient space RF /NΦ ); only when Σ ≻ 0 (i.e. under Assumption 4.3) does V = RF . The condition-number upper bound + decreases strictly only if the envelope ratio ML+1 /m+ L+1 of a newly added marginal is smaller than the current upper bound; in general one can only guarantee that the identification rank is unchanged. Assumption 4.2 (Positive and nondegenerate marginals). Each marginal pair (al , bl ) satisfies: (i) all entries are strictly positive, al,i > 0, bl,j > 0 for all i, j; (ii) both sides have full support (at least K distinct source and target states). Under these conditions the transport polytope M(al , bl ) has full tangent dimension (K − 1)2 (the generic case; the degenerate exceptions are point masses or zero-mass states, which are excluded by strict positivity). This is the condition under which the kernel of the Sinkhorn map’s Jacobian equals the gauge subspace G: ker J (l) = {∆θ : Φ(∆θ) ∈ G}, so a single marginal pair already reveals the full off-gauge 12

sensitivity. In practice condition (i) is verifiable by inspecting the data marginals; condition (ii) holds automatically for any nontrivial distribution. Randomly generated marginals (as in the experiments of this paper) satisfy this assumption with probability 1. F Assumption 4.3 (Linear independence off the gauge). The feature vectors {φ⊥ k }k=1 are linearly independent; equivalently rank(Σ) = F .

Proposition 4.4 (Global injectivity and the marginal lower bound). the positive-marginal P Under (l)⊤ J (l) is independent J assumption, the rank of the stacked sensitivity Gram matrix SL = L l=1 of the number of marginals L: for each l, ker J (l) = {∆θ : Φ(∆θ) ∈ G}, so rank(SL ) = rank(Σ) for every L ≥ 1. Identifiability requires only L ≥ 1; adding marginals does not change the rank. Identifiability necessarily requires that the off-gauge feature dimension not exceed (K − 1)2 . The role of the marginals is well-posedness, not identifiability. Remark 4.5 (The spectral scale of the Sinkhorn linearization).pThe core spectral bound (3.3) of Section 3 yields the lower bound σmin (J ) ≥ (πmin /(amax ε)) λmin (Σ) > 0 on the off-gauge subspace (derivation chain: Hessian ε diag(1/π) → inverse eigenvalue lower bound πmin /ε → composition with the feature Gram). When the unregularized optimal plan contains zero entries, πmin (ε)/ε may degenerate, causing the spectral lower bound to degenerate; monotonicity of the actual sensitivity Gram requires additional structural conditions on the cost, the marginals, and the features (numerical verification in Section 9).

4.2

Numerical verification

See Section 9 for details.

5

T2: Sparsistency

5.1

Statement

Setting. The true parameter θ ∗ ∈ RF is sparse, with support S = {k : θk∗ ̸= 0}. Write θbnlasso for the ℓ1 -penalized cross-entropy estimate, from which the support estimate Sbn is obtained by a relative threshold; an unpenalized refit on Sbn then gives the debiased estimate θeS,n . Theorem 5.1 (Sparsistency). Under the identifiability condition (Assumption A1) and an irrepresentability-type condition (Assumption A2), the support estimate Sbn and the debiased estimate θeS,n satisfy P(Sbn = S) → 1 (n → ∞), with exponentially decaying failure probability. Take a tuning sequence λn → 0, tn = O(λn ), nt2n → ∞. For a single marginal pair, or a single normalized empirical average, the scoreconcentration part of the failure probability is P(Sbn ̸= S) ≤ C1 exp(−C2 nt2n ),

C2 = 2/∆2max ,

where ∆max := maxk=1,...,F ∆k is the largest range of the per-sample score over all coordinates (the all-coordinate event (14) controls the active and the inactive coordinates simultaneously). At a fixed t one obtains a fixed exponential rate, but then λn cannot tend to zero; consistency 2 requires t = tn → 0, and the exponential rate becomes e−C2 ntn . If there are multiple independent marginal blocks with weights wl , the effective concentration scale of (15) should be used instead of n and a single ∆max directly. The Bernstein refinement is ! nt2n exp − 2 2 . 2σ + 3 ∆max tn 13

Accordingly, Sbn recovers S consistently, and the unpenalized refit θeS,n under correct specification is asymptotically unbiased on S. If one studies ε′ ̸= ε, the fixed-scale constraint of T3 and a support-stability assumption are additionally needed. Score range and the concentration inequality. The per-sample score sk = −∂θk log πij has mean zero under correct specification. The spectral sandwich (10) of the Sinkhorn linearization yields the absolute-value envelope for the range ∆k := max sk − min sk : ∆max ≤ 2Mmax = 2

πmax πmax max ∥PT φk ∥2 ≤ 2 max ∥φk ∥2 , επmin k=1,...,F επmin k=1,...,F

where the factor 2 comes from the two-sided values of the score (sk ∈ [−Mk , Mk ]); the range must not be underestimated as half the envelope. This gives the exponential rate constant C2 = 2/∆2max (with t2n kept in the exponent, free of the threshold). The Bernstein refinement C2B = 1/(2σ 2 + 23 ∆max tn ) is available when the score variance is estimable, but numerically the range bound at a fixed threshold t is more conservative (6×–156×), and the Bernstein variancedependent rate is not attained in practice. ∆max is a posteriorly estimable constant: once the plan is fitted, πmin , πmax , and the tangent-space projection can be computed to give a numerical estimate. Fine-grained concentration: conditioning on row counts. The analysis above treats all rows with the worst-case ∆max , which is a conservative upper bound. A more refined analysis can condition on the source-state counts N = (N1 , . . . , NK ): given N , the observations within each row are i.i.d., and the per-sample score ranges ∆i,k may differ across rows. The weighted Hoeffding bound conditional on N gives !
2n2 t2n ∗ ∗ P |∇ℓn,k (θ ) − ∇ℓk (θ )| ≥ tn | N ≤ 2 exp − PK , 2 i=1 Ni ∆i,k where ∆i,k is the range of the P per-sample score sk in row i. When Ni is close to its expectation nai , the denominator ≈ n i ai ∆2i,k , and the exponential rate is ∝ nt2n . A Chernoff bound controls the probability that some row count is too small: K  na 
X i P ∃i : Ni ≤ 12 nai ≤ . exp − 8 i=1

Under the total-probability decomposition P(failure) ≤ P(failure | N normal) + P(N abnormal), the Chernoff term affects only the leading constant C1 , not the exponential rate C2 . This refined analysis yields tighter constants when the score ranges differ substantially across rows, while the asymptotic order nt2n remains unchanged. Conditional assumptions (sufficient conditions for T2). The following are the standard sufficient conditions for Lasso support recovery [30], [26]; we state them uniformly on the OT information matrix H = ∇2 ℓ(θ∗ ) (the same matrix that connects to the population curvature of T4). A2.1 (Local curvature/RSC). There exist a neighborhood U ⊂ Θ of the true parameter and κ > 0 such that v ⊤ ∇2 ℓ(θ)v ≥ κ∥v∥2 for θ ∈ U , and the empirical Hessian satisfies a standard local RSC condition on this neighborhood. An additional requirement: the neighborhood is large enough that both the candidate solution constructed by the primal-dual witness and the global solution of the ℓ1 -penalized estimator lie in U ; if this fails, the conclusion below only guarantees the existence of a local solution satisfying the KKT conditions, not the global ℓ1 estimator. If the theorem’s conclusion is to concern the global ℓ1 estimator, a global selection condition must be added: for example, the penalized objective is convex on the entire parameter domain, or the objective values outside U are strictly higher than at the PDW candidate inside U . Local RSC alone does not exclude a lower global minimum outside the domain. 14

−1 A2.2 (Irrepresentability). For the support block decomposition of H = ∇2 ℓ(θ∗ ), ∥HS c S HSS ∥∞ ≤ 1 − η (η ∈ (0, 1)). −1 A2.3 (Beta-min and the penalty window). Let AS := HSS , cS := 2∥AS ∥∞ , and require e e SS − HSS ∥∞ ≤ δS the integrated matrix HSS of the restricted empirical Hessian to satisfy ∥H −1 1 e with ∥AS ∥∞ δS ≤ 2 . The Neumann series then gives ∥HSS ∥∞ ≤ ∥AS ∥∞ /(1 − ∥AS ∥∞ δS ) ≤ cS . The beta-min condition at the truth, mink∈S |θk∗ | > dS,n , and the dual-margin condition (2 − η)tn + rD,n < ηλn guarantee support and sign recovery, where dS,n := cS (λn + tn ) and rD,n is the local Taylor remainder of the inactive gradient. The primal-dual witness argument. On the all-coordinate event n o (14) Etn := ∇ℓn (θ∗ ) − ∇ℓ(θ∗ ) ∞ ≤ tn

(which controls the active and inactive coordinates simultaneously; controlling only the inactive coordinates does not control the active error, because the KKT expansion on the active coordinates contains the term ∇S ℓn (θ∗ ) − ∇S ℓ(θ∗ ), which is not constrained by an inactive-only event). First restrict the problem to the support coordinates S and set the inactive coordinates e + λn zS = 0 (zS ∈ ∂∥θeS ∥1 ). The to zero; the KKT equation on the active coordinates is ∇S ℓn (θ) ∗ e SS , together with the inTaylor expansion around θ , expressed with the integrated Hessian H vertibility of A2.3 and the all-coordinate event, gives ∥θeS −θS∗ ∥∞ ≤ cS (λn +tn ) = dS,n ; beta-min guarantees that the active coordinates do not cross zero. For the inactive coordinates, the KKT e < λn ; its leading term is controlled by HS c S H −1 , irrepresentability requirement is |∇S c ℓn (θ)| SS (A2.2) leaves a dual margin of at least η, and the local remainder is controlled by rD,n . Hence the constructed restricted candidate satisfies the full KKT conditions, giving support and sign recovery. Note: this PDW construction only guarantees that the candidate is a local solution satisfying the KKT conditions; for the candidate to equal the global ℓ1 -penalized estimator, one additionally assumes in A2.1 that the neighborhood U contains the PDW candidate and that the objective is sufficiently convex within U to exclude other global minima. Concentration of the all-coordinate score follows from Hoeffding’s inequality with a union bound over all F coordinates, giving the failure probability P(Etcn ) ≤ 2F exp(−2nt2n /∆2max ), where n is the number of independent samples per marginal pair. Once the support is obtained, the unpenalized refit on √ d −1 −1 S gives the debiased estimate, which satisfies n(θeS,n − θS∗ ) − ), where VSS is VSS HSS → N (0, HSS the score covariance. Here “debiased” means removing the shrinkage bias of the ℓ1 penalty, not strict finite-sample unbiasedness; a nonlinear M -estimator is generally only asymptotically unbiased. If the fitted ε′ differs from the generating ε, the residual entropic bias biasε = O(|ε′ − ε|) is given by T3 and does not vanish with n. Note also that under misspecification (ε′ ̸= ε) the support of the pseudo-true parameter θ † (ε′ ) may differ from S; an O(|ε′ − ε|) parameter deviation does not automatically imply support stability, which requires the additional assumption that beta-min exceeds the misspecification deviation. The concentration scale with multiple weighted marginals. P If the empirical gradient is a weighted sum of independent marginal blocks, Gn,k − Gk = L l=1 wl (Z l,k − EZl,k ), where the l-th marginal has nl samples and per-sample score range ∆l,k , then the concentration scale of the per-sample score contains at least L X wl2 ∆2l,k l=1

nl

,

(15)

one cannot simply take the smallest sample size without accounting for the marginal weights, nor use a single ∆max without marginal indices as the exact rate constant. The main setting of this paper is a single marginal pair (L = 1) or an equally-weighted pooling, in which case the above reduces to ∆2k /n; the general multi-marginal case requires rescaling the rate constant according to (15). 15

Lemma 5.2 (Mutual incoherence implies Gram-matrix irrepresentability). Let φk := vec(Φk ) and φ⊥ k := PT φk be the gauge-cleaned feature vectors. Define the mutual incoherence between the active set S and the inactive set S c as µ :=

max

k∈S c , j∈S

⊥ ⟨φ⊥ k , φj ⟩ .

Then the irrepresentability condition for the Gram matrix Σ is ΣS c S Σ−1 SS ∞ ≤

s3/2 µ < 1 − η, κS

κS := λmin (ΣSS ),

which holds when µ < (1 − η) κS /s3/2 (with s = |S|, η ∈ (0, 1)). This condition is a priorstructural sufficient condition, depending only on the feature matrices Φk and the gauge projector PG ⊥ , both known before estimation. Note that the conclusion of this lemma is limited to the irrepresentability of the feature Gram matrix Σ; the matrix used in Assumption A2.2 of T2 is H = ∇2 ℓ(θ ∗ ) (the OT information matrix), which differs from Σ under a general parameterization. The IR condition on Σ transfers to H only when there is a proportionality relation, a block-structure relation, or a direct matrix-comparison theorem between H and Σ; this bridge is not separately assumed in this paper. Hence the lemma gives a sufficient condition verifiable a priori from the feature structure, while the actual A2.2 still needs to be diagnosed post-fitting through an estimate of H.

5.2

Numerical verification

See Section 9 for details.

6

T3: Well-Posedness and Stability

6.1

Statement

Theorem 6.1 (Stability of the inverse map: local and global). Let Ψ be the inverse map from the observed conditional transition operator to the cost parameter, restricted to the identifiable quotient space RF /NΦ . Under identifiability (Assumption 4.3) and the positive-marginal assumption (Assumption 4.2), two statements hold. (a) Local strong monotonicity (no compactness needed). The feature-moment map M (θ) := Φ⊤ xθ ∈ RF , where xθ = vec(πθ ), is strongly monotone on every compact convex subset K ⊂ RF /NΦ (πmin (K) > 0, λmin (Σ) > 0): for any θ, θ ′ ∈ K, ⟨M (θ) − M (θ ′ ), θ − θ ′ ⟩ ≥ γK ∥θ − θ ′ ∥2 ,

γK := inf σmin (∇Mθ ) ≥ θ∈K

πmin (K) λmin (Σ) > 0. ε

The Jacobian ∇Mθ = Φ⊤ ∂x/∂θ has dimension F × F , matching the parameter dimension. −1 The moment-map inverse ΨM is γK -Lipschitz on the image of K; composing with q 7→ x = Sa q ⊤ and Φ yields the Lipschitz bound of the Q-space inverse map Ψ (see the global part below). (b) Global strong monotonicity (compact domain required). To obtain a uniform Lipschitz constant over the entire parameter domain, assume Assumption 6.6 (Θ ⊂ RF /NΦ is compact and convex, πmin := inf θ∈Θ minij πij (θ) > 0). Σ = Φ⊤ PT Φ is a constant matrix, and λmin (Σ) > 0 is guaranteed by T1. Then the bound in (a) holds uniformly over Θ: ε πmin λmin (Σ) > 0, L := γ −1 ≤ , γ := inf σmin (∇Mθ ) ≥ θ∈Θ ε πmin λmin (Σ) The moment-map inverse ΨM is L = γ −1 -Lipschitz. Composing xθ = Sa qθ (Sa = diag(a1 , . . . , aK )⊗ IK ) with the feature reweighting Φ⊤ Sa , the Lipschitz bound of the Q-space inverse map Ψ is LΘ := ∥Φ⊤ Sa ∥op γ −1 ≤ 16

ε ∥Φ⊤ Sa ∥op . πmin λmin (Σ)

The promotion from pointwise to global relies on the compactness of Θ (the conditioning hierarchy of Section 3); without compactness there is only the pointwise lower bound and strong monotonicity degrades to local strong monotonicity. The unbounded case (global strong monotonicity on RF /NΦ ) is an open problem. In practice Θ can be enforced by projected gradient descent or by a prior; the algorithm need not know the boundary of Θ to apply the local result during optimization. Theorem 6.2 (Entropic bias). When the data are generated at level ε and the estimator is fitted at level ε′ , work at the population level (infinite-sample limit) under a fixed smooth scale constraint g(θ) = 0 (e.g. g(θ) = (∥θ∥2 − 1)/2). Define the pseudo-true branch θ † (ε′ ) and its scale multiplier ν † (ε′ ) as the unique local branch of the augmented estimating equation   ∇θ CE(Qθ∗ ,ε ∥Qθ,ε′ ) + ν∇g(θ) F (θ, ν, ε′ ) := =0 g(θ) Here “unique” requires the augmented Jacobian to be invertible, which is an implicit function theorem condition. Assume this branch exists and that the augmented Jacobian D(θ,ν) F is invertible and continuous along the branch. By the implicit function theorem, d(θ † , ν † ) = −[D(θ,ν) F ]−1 ∂ε′ F . dε′ Let q † (ε′ ) := Qθ† (ε′ ),ε′ ,

q ∗ := Qθ∗ ,ε .

Then the total derivative of the conditional operator along the branch is dθ † dq † = J + ∂ε′ Q θ† ,ε′ , † θ dε′ dε′ where the direct regularization term is given by Proposition 6.3 and q = Ra x. If this total derivative is bounded in a neighborhood of ε, then biasε′ := ∥q † (ε′ ) − q ∗ ∥F ≤ c |ε′ − ε|, and the bias vanishes at ε′ = ε. The Jacobian of the augmented estimating equation is a constrained version of the cross-entropy Hessian, distinct from the Jacobian of the moment map M . Proposition 6.3 (Differentiability of the Sinkhorn plan in ε). Fix marginals (a, b) and a cost C. The entropic OT plan πε = πε (C, a, b) is differentiable in ε > 0. Write xε = vec(πε ), let 2 2 B ∈ RK ×(K−1) be an orthonormal basis of the tangent space T , and HT = εB ⊤ Dπ−1 B the restricted Hessian (Lemma 3.1). Then the exact derivative is ∂ε xε = −BHT−1 B ⊤ log xε , satisfying the norm bound (the spectral sandwich (10) guarantees that the exact form and the spectral proxy share the same upper and lower bounds) ∥∂ε πε ∥F ≤

πmax πmax ∥PT (log πε )∥F ≤ K log(1/πmin ), ε ε

where πmax := maxij πij (ε), πmin := minij πij (ε). Proposition 6.4 (Local minimum of the entropic bias function). Under the fixed-scale and augmented implicit-function regularity conditions of Theorem 6.2, define B2 (ε′ ) := ∥q † (ε′ ) − q ∗ ∥2F . 17

If q † is twice differentiable in a neighborhood of ε, then B2′ (ε) = 0 and 2

B2′′ (ε) = 2

dq † (ε) ≥ 0. dε′ F

If one further assumes dq † /dε′ (ε) ̸= 0, then B2′′ (ε) > 0, so that ε′ = ε is a strict local minimum of the squared bias, and hence a strict local minimum of the original bias. The original Frobeniusnorm bias is typically not differentiable at zero, which is why the squared bias is used. Remark 6.5 (Empirical comparison). The measured local inverse-map Lipschitz constant is Lemp = 25.8. The theoretical prediction LΘ = ε∥Φ⊤ Sa ∥op /(πmin λmin (Σ)) is derived by composing the core spectral bound (3.3) with the feature reweighting; using the πmin of the fitted plan and ∥Φ⊤ Sa ∥op gives a data-dependent prediction. Derivation. ∇Mθ = Φ⊤ Sπθ Φ = Φ⊤ BHT−1 B ⊤ Φ is symmetric positive semidefinite. By the spectral sandwich (10), for any v ∈ RF , v ⊤ ∇Mθ v = (B ⊤ Φv)⊤ HT−1 (B ⊤ Φv) ≥

πmin (θ) πmin (θ) ∥PT Φv∥2 ≥ λmin (Σ)∥v∥2 . ε ε

Note that the strong-monotonicity modulus carries λminp (Σ) to the first power, whereas the singular-value bound (3.3) of the Jacobian Jθ carries λmin (Σ); the two scales differ. On a convex domain (K or Θ), the segment tθ + (1 − t)θ ′ stays in the domain, so the uniform bound applies along the entire segment; the fundamental theorem of calculus gives strong monotonicity, and Cauchy–Schwarz gives the reverse Lipschitz property of the moment map, so ΨM is γ −1 -Lipschitz. The Q-space inverse map Ψ is obtained through q 7→ x = Sa q (Sa = diag(a1 , . . . , aK ) ⊗ IK = Ra−1 ) composed with Φ⊤ and then ΨM . The monotone-operator viewpoint. On the quotient space, M : θ 7→ Φ⊤ xθ is a strongly monotone operator. By the Zarantonello–Minty theorem [29, 17], a strongly monotone map is injective with a Lipschitz inverse; [21] gives the proximal-point algorithmic framework for monotone operators. The paper-specific content is the uniform lower bound on the modulus γ derived from the Sinkhorn linearization (Section 3). Assumption 6.6 (Compact domain for the global statement). This assumption is needed only for the global statement T3(b); the local statement T3(a) does not require it. Assume the parameter is restricted to a compact convex domain Θ ⊂ RF /NΦ on which the Sinkhorn plan entries are uniformly bounded away from 0: πmin := inf min πij (θ) > 0, θ∈Θ ij

and the gauge-cleaned Gram has a minimal eigenvalue: λΘ min := λmin (Σ) > 0 (the global version of the pointwise identifiability condition).

6.2

Numerical verification

See Section 9 for details.

7

T4: Convergence

7.1

Statement

Theorem 7.1 (Monotone convergence to a local minimum). Under correct specification, let the population cross-entropy objective be ℓ(θ) =

L X

(l)

wl CE(Qθ∗ ∥Qθ (al , bl )),

l=1

X l

18

wl = 1, wl ≥ 0,

(l)

where Qθ∗ is the conditional transition operator at the true parameter θ ∗ (the population objective, not the empirical objective ℓn ). The population Hessian at the true parameter θ ∗ is positive definite (derived from the core spectral bound), and gradient descent with a sufficiently small step size converges monotonically to a local minimum, provided all iterates remain within the local strong-convexity neighborhood of θ ∗ (e.g. the sublevel set containing the initial point). For the statistical theory of entropic OT see [5] (central limit theorems), [16] (statistical bounds), and [20] (entropy estimation).

7.2

Proof roadmap

Population cross-entropy andPweights. To cover multiple marginals simultaneously, let l = 1, . . . , L with weights wl ≥ 0, l wl = 1; for each source state i let the cross-entropy weight be ρl,i > 0. The operator-level objective of this paper weights every row equally, i.e. ρl,i = 1; a likelihood weighted by the joint sample frequencies corresponds to general ρl,i . The population loss is L K K X X X (l) (l) ℓ(θ) = − wl ρl,i Qθ∗ (i, j) log Qθ (i, j). l=1

i=1

j=1

Exact form of the Hessian (the second-order residual vanishes). Write Jl (θ) = (l) (l) ∂qθ /∂θ and the diagonal weight at the truth Wl = diagi,j (ρl,i /Qθ∗ (i, j)). Differentiating the P population loss once gives ∇ℓ(θ) = − l wl Jl (θ)⊤ (ρl ⊙ ql∗ /ql (θ)) (division coordinatewise). Differentiating again, two types of terms appear at θ ∗ : the Jacobian outer product, and the second-derivative term of ql . For a fixed source state i, the coefficient of the second-order P (l) P (l) term is −wl ρl,i j ∂ 2 Qθ (i, j)/∂θ∂θ ⊤ θ∗ . Since j Qθ (i, j) = 1 holds identically in θ, double P (l) differentiation gives j ∂ 2 Qθ (i, j)/∂θ∂θ ⊤ = 0, so the second-derivative residual term vanishes row by row, with no cross terms. The population Hessian at the truth is exactly the Fishertype matrix L X 2 ∗ ∇ ℓ(θ ) = wl Jl (θ ∗ )⊤ Wl Jl (θ ∗ ) . l=1 (l)

Curvature lower bound. From Qθ∗ (i, j) ≤ 1 we get Wl ⪰ ρmin I, where ρmin := minl,i ρl,i . Combined with the core spectral bound (3.3), λmin (∇2 ℓ(θ ∗ )) ≥

X

wl ρmin σmin (Jl (θ ∗ ))2 ≥

X

l

wl ρmin



 πmin,l 2 λmin (Σ). al,max ε

l

For the operator-level objective, take ρmin = 1 and use the looser bound amax ≤ 1, giving λmin (∇2 ℓ(θ ∗ )) ≥

2 πmin,∗

ε2

λmin (Σ) =: µ > 0,

πmin,∗ := min πmin,l . l

When the right-hand side is strictly positive, ℓ has a positive definite Hessian at the truth, derived from the core spectral bound. By continuity of the Hessian, there exist a neighborhood U of the truth and constants 0 < µ0 < µ and Lℓ < ∞ such that µ0 I ⪯ ∇2 ℓ(θ) ⪯ Lℓ I (θ ∈ U ), so ℓ is locally strongly convex on U with Lipschitz gradient. µ0 is generally strictly smaller than µ (by Hessian continuity, µ0 = µ − δ for some δ > 0). Gradient descent θr+1 = θr − η∇ℓ(θr ), when the initial point lies in the basin of attraction within U and 0 < η ≤ 1/Lℓ , satisfies ℓ(θr+1 ) − ℓ(θ ∗ ) ≤ (1 − ηµ0 )(ℓ(θr ) − ℓ(θ ∗ )), converging monotonically to the local minimum. Proposition 7.2 (Concentration of the empirical Hessian). Let θbnCE be the empirical minimizer of the normalized empirical cross-entropy objective, distinct from the ℓ1 -penalized estimate θbnlasso

19

and the debiased refit estimate θeS,n of T2. On a neighborhood U of the true parameter, assume the full per-sample cross-entropy Hessian contribution ∇2θ ψ(Z; θ) satisfies sup ∥∇2θ ψ(z; θ)∥op ≤ MH ,

θ∈U,z

and that the empirical Hessian is Lipschitz in the parameter: ∥∇2 ℓn (θ) − ∇2 ℓn (θ ′ )∥op ≤ LH ∥θ − θ ′ ∥. If for some rn > 0 one has the high-probability estimator error bound   P ∥θbnCE − θ ∗ ∥ > rn ≤ δθ (n), then for any t > 0, matrix Hoeffding and the triangle inequality give     nt2 2 CE 2 ∗ b P ∥∇ ℓn (θn ) − ∇ ℓ(θ )∥op > t + LH rn ≤ 2F exp − + δθ (n). 2 2MH In particular, if LH rn ≤ µ/4 and one takes t = µ/4, then whenever n≥

2 32MH 2F log , µ2 δH

δH + δθ (n) ≤ δ,

with probability at least 1 − δ one has   µ λmin ∇2 ℓn (θbnCE ) ≥ > 0. 2 2 Here µ = (πmin,∗ /ε2 )λmin (Σ) is the curvature lower bound of the population Hessian. The rate rn = O(n−1/2 ) requires a separate proof of the corresponding estimator error tail bound.

The IOT objective has no known analogue of a condition that rules out spurious local minima; this paper retains the local-minimum structure.

7.3

Numerical verification

See Section 9 for details.

8

O5: Misspecification Analysis

8.1

Statement

Observation 8.1 (Projection onto the OT model set). Let MOT be the model set of conditional transition operators generated by entropic OT, and Q⋆ the true (possibly non-OT) transition operator. Under a compact convex parameter domain Θ (Assumption 6.6), the image of Qθ is compact; define the projection residual Res(Q⋆ ) = minQ′ ∈MOT CE(Q⋆ ∥Q′ ) (the minimum is attained). Without compactness, replace min by inf. (1) Convergence (standard M -estimation). The IOT estimator converges to the pre-image of the projection of Q⋆ onto MOT . If the projection minimizer is unique and lies in the identifiable region R, then   p θbn → − Ψ arg ′min CE(Q⋆ ∥Q′ ) , Q ∈MOT

where Ψ is the inverse map of T3 (well-defined and Lipschitz on the identifiable region R by T1 and T3). If the projection minimizer is not unique (multiple minima), the statement can only be phrased as convergence in probability of θbn to the set of minimizers in distance; if the 20

minimizer falls outside R, then Ψ is set-valued (gauge-kernel ambiguity) and the conclusion is set-valued convergence. (2) Residual bound (definitional). The limiting residual Res(Q⋆ ) is nonnegative, and for any random OT model Qrand ∈ MOT lying on the manifold, Res(Q⋆ ) ≤ CE(Q⋆ ∥Qrand ), a bound that is the definitional minimality of the projection: Res does not exceed the value at any feasible point. Remark 8.2 (Dependence of O5 on T1–T4 and the Hölder exponent). O5 depends on T1 (the identifiable region) and T3 (the Lipschitz inverse map), and not on T2 or T4. The Hölder continuity of the projection map Q 7→ arg minQ′ ∈MOT CE(Q∥Q′ ) controls the robustness of the estimator under misspecification. The experiments (Section 9) estimate the empirical effective exponent αeff ∈ (0, 1) via log-log regression across multiple settings; the ε = 0.1 setting is affected by Adam optimization residual, and the reliability of αeff should be assessed jointly with the condition number analysis (T1) and optimization convergence (T4).

8.2

Proof roadmap (consistency)

The pseudo-true projection. Let the true conditional transition operator be Q⋆ (not necessarily from the OT model of this paper), the model manifold MP {Qθ : θ ∈ Θ}, and the OT = P population conditional cross-entropy LQ⋆ (θ) = CE(Q⋆ ∥Qθ ) = − i ρi j Q⋆ (i, j) log Qθ (i, j) (ρi > 0 the source-state weights). The operator-level objective of this paper takes ρi = 1 (equal row weights), consistent with ρl,i = 1 in the T4 population objective above; if the samples are drawn at the joint frequencies, ρi can be taken as the empirical marginal frequency ai of the source states. If Q⋆ ∈ MOT and T1 holds, the projection is Q† = Q⋆ ; otherwise Q† is the pseudo-true model point. The consistency argument. First make explicit the empirical process to which the emb is the row-normalized conditional operator and pirical objective corresponds. In this paper Q the empirical objective uses equal row weights i = 1), so the ULLN is stated for this operatorP P (ρ b ⋆ level empirical process Ln,Q (θ) := − i j Qn (i, j) log Qθ (i, j); if the samples are drawn at the joint frequencies and weighted by ρi = ai , the empirical process is different, and the corresponding weighted ULLN must be proved separately—the equal-row-weight form here cannot be applied directly. Assume Θ is compact, Qθ (i, j) is continuous and uniformly positive on Θ, and the population objective has a unique minimizer θ † ∈ arg minθ∈Θ LQ⋆ (θ). Then the function family {− log Qθ (i, j) : θ ∈ Θ} is uniformly bounded and equicontinuous on the finite state space, and the operator-level empirical objective Ln,Q⋆ satisfies the uniform law of large numbers p sup Ln,Q⋆ (θ) − LQ⋆ (θ) → − 0. θ∈Θ

Take any convergent subsequence θbnr ; its limit θ̄ satisfies LQ⋆ (θ̄) ≤ LQ⋆ (θ) (θ ∈ Θ) by empirical minimality and the ULLN; unique minimization gives θ̄ = θ † . All convergent subsep p quences have the same limit, so θbn → − θ † and Qθbn → − Q† . If the projection is not unique (the population objective has multiple global minima), the statement can only be phrased as convergence of the distance from the estimator to the set of minimizers S † := arg minθ∈Θ LQ⋆ (θ): p dist(θbn , S † ) → − 0. If θbn and θ † both lie in the same compact convex regular domain Θreg and Q† lies in its model image, T3 gives the local stability transfer ∥θbn − θ † ∥ ≤ LΘ ∥Qθbn − Q† ∥.

8.3

Numerical verification

See Section 9 for details.

21

9

Experiments

Protocol. All experiments fix the numpy and PyTorch random seeds for reproducibility; the main setting is ε = 0.2, K ∈ {4, 5, 6}, and the number of features F ∈ {6, 8, 10, 15}, with E6/E7/E10 additionally scanning ε; marginals are generated randomly and the sample size ranges from 1000 to 40000; Sinkhorn plans are computed in double precision. The selfcontained implementation includes a differentiable Sinkhorn, the feature-parameterized cost, Adam with cosine annealing, gradient clipping, and non-finite restart. The aggregated results are in 09 iot theory/output/experiments e1 e10.json.

9.1

T1: Identifiability (E6, E7)

Setting Full rank Collinear (feature 2 ≈ feature 1) Pure gauge feature

Rank

Condition number

Direction recovery correlation

6 6 5

2.58 2474 —

0.999 0.958 not recoverable

Table 2: Comparison of identifiability settings Verification of the marginal lower bound: K=4, F =6 already reaches full rank 6 at L = 1; K=3, F =5 is capped at (K − 1)2 = 4 and K=4, F =15 is capped at 8, and adding L has no effect; for K=5, F =10 the condition number improves with L from numerical ill-posedness (with a single marginal S1 is nearly singular, with condition number ≈ 1016 , recorded as ∞) down to 2582—an improvement in numerical well-posedness, not a repair of the rank (T1 has proved that the rank is independent of L). Collinearity damages the condition number and unique identifiability (cond 2.58 → 2474), while the direction remains recoverable (corr 0.958); the pure gauge feature is the only true failure point. Condition-number divergence as ε → 0 (E6, 10 seeds). For ε ∈ {0.5, 0.3, 0.2, 0.1, 0.05, 0.02} the log means are log cond(SL ) ∈ {5.19, 6.83, 3.53, 4.82, 10.77, 15.75}, rising by about 10.6 log units (a factor of about 3.8 × 104 ) from ε = 0.5 to ε = 0.02; the small-ε tail shows a divergent trend for the given features and marginals. The middle segment (0.2, 0.1) fluctuates because random features occasionally come close to the gauge directions, showing that the actual condition number need not be strictly monotone. This corresponds to the degeneration of πmin → 0 in the core spectral bound (3.3); the trend is an empirical observation for the given settings.

Generality (E7). The divergence trend is reproduced under (K, F, ε) ∈ {(6, 8, 0.2), (4, 6, 0.2), (8, 12, 0.2), (6, 8 (from ε = 0.5 to 0.05, cond rises 699 → 5.8 × 107 , 3.3 × 104 → 1.2 × 105 , 478 → 9.5 × 104 , 444 → 5.4 × 104 , and 342 → 1.3 × 107 respectively).

22

Figure 1: Condition-number divergence as ε → 0

Figure 2: Generality verification

23

9.2

T2: Sparsistency (E1, E3, E8, E9)

Recovery curve (E1). The support recovery probability is estimated from 100 seeds per sample size over n ∈ {200, 500, 800, 1200, 2000, 3000, 5000, 8000, 12000, 20000}, with binomial 95% confidence intervals. Figure 3 shows a four-panel diagnostic: (a) sample-complexity phase transition, with recovery probability rising from 0.08 at n = 200 to 0.93 at n = 8000, with a slight dip to 0.90 at n = 12000 before reaching 0.96 at n = 20000; (b) concentration mechanism decomposition, where the score-event probability ∥Gn ∥∞ ≤ tn rises from 0.19 to 1.00 and the row-count event probability from 0.89 to 1.00; (c) failure-probability rate diagnostic, with the linear fit of − log Pbfail against nt2n /Amax verifying the exponential decay rate C2 of Theorem 5.1; (d) error decomposition, where the false-positive rate drops from 0.26 to 0.008 and the falsenegative rate from 0.05 to 0, confirming that support failure is primarily driven by false positives.

Figure 3: Theorem-compatible sparsistency diagnostics: (a) sample-complexity phase transition; (b) concentration mechanism decomposition; (c) failure-probability rate diagnostic; (d) decomposition of support failure causes

Voting (E3, 100 seeds). Each restart is run on a different bootstrap subsample to decorrelate the failures; with vote counts R ∈ {1, 3, 5, 7, 9} the recovery probabilities are {0.20, 0.40, 0.30, 0.60, 0.50}, and the single-restart success rate is 0.47. Voting raises the recovery rate from 0.47 to 0.50– 0.60 with fluctuations; voting is a conditional stabilizer, the gain being limited by a residual systematic failure component. Numerical verification of irrepresentability (E8b, 40 seeds). The irrepresentability quantity is computed on the population OT information matrix H = ∇2 ℓ(θ∗ ) (averaged over L = 8 marginal pairs), landing in [1.36, 7.43] with mean 3.19; the condition < 1 is satisfied 0/40 times. Although IR is not satisfied, recovery still succeeds in E1b, so the condition is positioned as a strictly sufficient one: it is generically violated on unstructured features, while the recovery conclusion of the theorem still holds empirically. Lemma 5.2 gives a sufficient condition verifiable 24

Figure 4: Effect of the voting mechanism on the recovery probability a priori from the feature structure, but that condition is stronger than directly checking H-IR; in the experimental setting of this paper, the sufficient condition of the lemma is not met, and the directly computed H-IR condition fails as well; consequently, the lemma does not apply in this setting. The failure of the sufficient condition while recovery still succeeds indicates that the condition, though sufficient, is not necessary in this setting.

Figure 5: Numerical verification of the irrepresentability condition

Comparison of Bernstein constants (E9, 8 seeds). The ratio of the Hoeffding constant C2 = 2/∆2max to the Bernstein refinement C2B lies in [13.7, 33.9] (mean 24.4): the variance-aware Bernstein constant is 1–2 orders of magnitude larger than the range-based Hoeffding constant, confirming that the refinement is substantially sharper.

25

Figure 6: Comparison of the Bernstein and Hoeffding constants

9.3

T3: Well-posedness (E2) Perturbation δQ 0.0003 0.001 0.003 0.01

Transfer error 0.0018 0.0030 0.0067 0.0076

Table 3: Perturbation transfer error The ε′ -bias V-shape (E2, 40 seeds). The transfer bias is averaged over 40 seeds with mean ± standard error reported; on the grid ε′ ∈ {0.05, 0.1, 0.2, 0.3, 0.5} it is 0.742 ± 0.037, 0.718 ± 0.046, 0.663 ± 0.020, 0.724 ± 0.032, 0.819 ± 0.047. The quadratic fit has significant positive curvature (p = 0.024, R2 = 0.77), so the V-shape is statistically supported; the empirical minimum is at ε′ ≈ 0.2, matching the generating ε = 0.2. In the experiments of this paper, the V-shape is an empirical observation that provides finite-sample support for the local minimality of the squared bias. The original Frobenius-norm bias B(ε′ ) is theoretically non-differentiable at the truth (the V-shaped cusp), and the V-shape observed in the experiment itself supports this fact; the squared bias removes the cusp issue.

26

Figure 7: The ε′ -bias V-shape

9.4

T4: Convergence (E4)

Comparison of initialization schemes (E4, 40 restarts). The success rates of the three schemes—random / multiscale (coarse warm-start) / prior-guided (initialized on the true support)— are 0.25, 0.95, 0.20 respectively, with mean cross-entropies 14.39, 14.34, 15.03. The multiscale coarse warm-start raises the success rate from 0.25 to 0.95 (about 3.8×) with the lowest objective: local strong convexity (Theorem 7.1) guarantees that fixed-step-size gradient descent decreases monotonically within a basin, and the multiscale scheme raises the fraction of initializations landing in the basin containing the global optimum. The actual optimization in this paper uses Adam, cosine annealing, and gradient clipping; their convergence is an empirical phenomenon—Theorem 7.1 addresses fixed-step-size gradient descent and does not cover adaptive optimizers such as Adam. The prior-guided scheme performs worst, because fixing the true support injects the correct active set but biases the magnitudes during optimization. Failures under random initialization are systematic: a restart stuck in a spurious basin does not escape by switching to a different random initialization.

27

Figure 8: Comparison of initialization schemes

9.5

O5: Misspecification (E5)

Under a non-OT generating mechanism (row softmax plus noise), the projection residual is 4.50, larger than the true-OT residual 2.43 but much smaller than the random-model residual 10.63 (ratio 0.42); the true-OT residual 2.43 is the total cross-entropy over all source states (about 0.4 nats per state), which under correct specification should approach the entropy H(Q∗ ) of the true distribution, with the gap to H(Q∗ ) reflecting finite-sample estimation error. The fitted θ collapses toward 0 (magnitude 0.03–0.09), indicating that the data cannot be explained by an OT cost. Hölder exponent across settings (E5). The Hölder exponent α of the projection map is estimated by a log–log regression of the fitted-parameter displacement against the observation perturbation, reproduced under four (K, F, ε) settings: αeff ∈ {0.353, 0.266, 0.369, 0.298} for (6, 8, 0.2), (5, 6, 0.2), (8, 12, 0.3), (6, 8, 0.1) respectively (per-setting R2 ∈ [0.53, 0.92], p ∈ [4 × 10−5 , 0.025]), all lying in (0, 1) and away from the degenerate endpoint 1. αeff is a settingdependent empirical exponent estimate; the cross-setting variation reflects the influence of ε and problem dimension on the condition number of the projection map. The ε = 0.1 setting yields αeff = 0.298, the lowest among the four. A quantitative analysis reveals that this setting is contaminated by the Adam optimization residual: comparing the projection results after 500 iterations (the experimental setting) and 2000 iterations (closer to the true minimum), the optimization residual is ∥θ500 − θ2000 ∥ = 0.267, which exceeds the true parameter displacement at the three smallest perturbation amplitudes (10−5 , 3 × 10−5 , 10−4 )— 0.107, 0.093, and 0.248, with residual ratios 2.51, 2.87, and 1.08, respectively. After removing these three residual-dominated points, the log-log fit over the remaining six signal points yields α = 0.048 (p = 0.67, R2 = 0.05), losing statistical significance. At ε = 0.1, the elevated condition number (T1) drives the projection map into an ill-conditioned regime; although the CE has flattened after 500 iterations (∆CE ≈ 10−5 ), the parameters continue to drift, and the optimization residual contaminates the log-log slope estimate at small perturbations. In the other three settings (ε ≥ 0.2), the optimization residual is only 0.003–0.016, and all perturbation responses are well above the residual (ratios < 0.42), so the α estimates are reliable.

28

Figure 9: Verification of the Hölder exponent across settings

Figure 10: Adam optimization residual analysis at ε = 0.1: the CE has flattened after 500 iterations while the parameters continue to drift; the responses at the three smallest perturbation amplitudes are submerged by the optimization residual

29

9.6

The ε-scaling law of πmin (E10)

The minimum plan entry πmin (ε) is measured at ε ∈ {0.05, 0.1, 0.2, 0.3, 0.5}, averaged over random marginals in log space. An empirical fit of πmin ≥ cεα over the finite window (log–log regression) gives the effective exponent αeff ≈ 19.9 (R2 = 0.89, p = 0.015). These exponents are empirical fits over a finite window, not global polynomial lower bounds; exponential decay e−∆/ε can also produce a large effective power exponent over a finite window, and the observed αeff ≈ 19.9 may well be exactly this phenomenon. In the unstructured case πmin collapses across many orders of magnitude (from 5.8 × 10−13 at ε = 0.5 to 2.3 × 10−31 at ε = 0.05). The large effective exponent αeff ≈ 19.9 places the random features in the severely degenerate regime. The associated Lipschitz constant L ≤ ε/(πmin p λmin (Σ)) inflates sharply in the degenerate regime; the effective information scale is πmin (θ, ε) λmin (Σ)/ε, and a threshold cannot be given in terms of cond(Σ)−1 alone. The statement ε ≫ cond(Σ)−1 in this paper is an empirical heuristic.

Figure 11: The ε-scaling law of πmin

10

Conclusion

This paper establishes the statistical and algorithmic theory of IOT under a feature-parameterized cost. The central technical contribution is the Sinkhorn linearization and the spectral proxy (Section 3): the spectral sandwich of the restricted Hessian yields a single core spectral bound, which, under the respective regularity conditions, drives Theorem 4.1 (identifiability), Theorem 6.1 (well-posedness), and Theorem 7.1 (population convergence with a conditional finitesample statement); Theorem 5.1 (sparsistency) additionally requires irrepresentability of the actual OT information matrix, all-coordinate score concentration, and a global selection condition; O5 (misspecification) interprets the estimation target as the pseudo-true projection onto the OT model set. The main open problems: minimax lower bounds for IOT estimation, global strong monotonicity without compactness, a proof of the Hölder-continuity conjecture, and extensions to continuous state spaces. All numerical verifications (E1–E10) are in Section 9.

References [1] F. Andrade, G. Peyré, and C. Poon, Sparsistency for inverse optimal transport, arXiv 30

preprint arXiv:2310.05461, (2024). [2] Z. Bao et al., Well-posedness of the inverse bregman optimal transport problem, (2025). Working paper. [3] D. Bartl, V. Chernozhukov, and J. Niles-Weed, Curvature of optimal transport with respect to the cost and applications to inverse optimal transport, arXiv preprint arXiv:2604.22670, (2026). [4] E. Bernton, P. Ghosal, and M. Nutz, Stability of entropic optimal transport and schrödinger bridges, Journal of Functional Analysis, 283 (2022), p. 109622. [5] J. Bigot, E. Cazelles, and N. Papadakis, Central limit theorems for entropyregularized optimal transport on finite spaces and statistical applications, Electronic Journal of Statistics, 13 (2019), pp. 5112–5150. [6] R. Cominetti and J. San Martı́n, Asymptotic analysis of the exponential penalty trajectory in linear programming, Mathematical Programming, 67 (1994), pp. 169–187. [7] M. Cuturi, Sinkhorn distances: Lightspeed computation of optimal transport, Advances in Neural Information Processing Systems, 26 (2013). [8] A. Dupuy and A. Galichon, Personality traits and the marriage market, Journal of Political Economy, 122 (2014), pp. 1271–1319. [9] A. Galichon, Optimal Transport Methods in Economics, Princeton University Press, 2018. [10] A. González-Sanz, M. Groppe, and A. Munk, Identifiability and exact reconstruction of the optimal transport cost on finite spaces, arXiv preprint arXiv:2410.23146, (2024). [11] A. González-Sanz, M. Groppe, and A. Munk, Nonlinear inverse optimal transport: identifiability of the transport cost from its marginals and optimal values, SIAM Journal on Mathematical Analysis, 56 (2024), pp. 7808–7829. arXiv:2312.05843. [12] M. Hallin, E. del Barrio, J. Cuesta-Albertos, and C. Matrán, Center-outward distribution and quantile functions, annular spreads, and the Mallows–Wasserstein distances between distributions, Annals of Statistics, 49 (2021), pp. 2567–2605. [13] A. Korba and A. Rudi, Statistical estimation of optimal transport maps: a density estimation approach to learning the optimal transport map between distributions, arXiv preprint arXiv:2106.07360, (2021). [14] C. Léonard, A survey of the schrödinger problem and some of its connections with optimal transport, Discrete & Continuous Dynamical Systems, 34 (2014), pp. 1533–1574. [15] M. Mascherpa, A. Ringh, A. Taghvaei, and J. Karlsson, A convex approach for Markov chain estimation from aggregate data via inverse optimal transport, arXiv preprint arXiv:2511.16458, (2025). [16] G. Mena and J. Niles-Weed, Statistical bounds for entropic optimal transport: sample complexity and the central limit theorem, in Advances in Neural Information Processing Systems, vol. 32, 2019. [17] G. J. Minty, Monotone (nonlinear) operators in Hilbert space, Duke Mathematical Journal, 29 (1962), pp. 341–346. [18] S. Persiianov et al., Semi-supervised inverse entropic optimal transport, (2025). Working paper. 31

[19] G. Peyré and M. Cuturi, Computational optimal transport, Foundations and Trends in Machine Learning, 11 (2019), pp. 355–607. [20] A.-A. Pooladian and J. Niles-Weed, Entropic estimation of optimal transport maps, arXiv preprint arXiv:2109.12004, (2021). [21] R. T. Rockafellar, Monotone operators and the proximal point algorithm, SIAM Journal on Control and Optimization, 14 (1976), pp. 877–898. [22] F. Santambrogio, Optimal Transport for Applied Mathematicians, vol. 87 of Progress in Nonlinear Differential Equations and their Applications, Birkhäuser, 2015. [23] R. Sinkhorn, A relationship between arbitrary positive matrices and doubly stochastic matrices, The Annals of Mathematical Statistics, 35 (1964), pp. 876–879. [24] A. M. Stuart and M.-T. Wolfram, Inverse optimal transport, SIAM Journal on Applied Mathematics, 80 (2020), pp. 599–619. [25] A. Vacher, B. Muzellec, A. Nishimura, and A. Dunipace, Bayesian inverse optimal transport, in Proceedings of Machine Learning Research (PMLR), 2021. arXiv:2105.12140. [26] M. J. Wainwright, Sharp thresholds for high-dimensional and noisy sparsity recovery using ℓ1 -constrained quadratic programming (Lasso), IEEE Transactions on Information Theory, 55 (2009), pp. 2183–2202. [27] J. Weed, An explicit analysis of the entropic penalty in linear programming, in Conference on Learning Theory (COLT), 2018. [28] H. White, Maximum likelihood estimation of misspecified models, Econometrica, 50 (1982), pp. 1–25. [29] E. H. Zarantonello, Solving functional equations by contractive averaging, Tech. Report Technical Report #160, Mathematics Research Center, U.S. Army, 1960. [30] P. Zhao and B. Yu, On model selection consistency of Lasso, Journal of Machine Learning Research, 7 (2006), pp. 2541–2563.

32

Related documents

• Recent Advances in Optimal Transport for Machine Learning#350625
• A Gauge-Fixed Sinkhorn Algorithm for Semi-Relaxed Optimal Transport#246916
• Variational objectives for amortized Bayesian inference in inverse problems: The role of posterior conditioning#1051638
• Inverse scattering in biological samples via beam propagation.#479318
• Sobolev Gradient Ascent for Optimal Transport: Barycenter Optimization and Convergence Analysis#153063