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

The paper develops the statistical and algorithmic theory of inverse optimal transport (IOT) under the feature-parameterized cost Cθ(i, j) = −θ⊤φ(i, j), where observations are conditional transition operators of the entropic optimal transport plan. The central 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 (preserves all singular-value bounds) yet geometrically transparent.

The restricted Hessian of entropic OT on the tangent space satisfies the spectral sandwich (πmin/ε)I ⪯ HT−1 ⪯ (πmax/ε)I, from which the single core spectral bound σmin(Jθ) ≥ (πmin/(amax ε))√λmin(Σ) is derived, driving the entire theory.

T1 (Identifiability): θ is identifiable on the quotient of the gauge kernel (RF/NΦ), with dimension bound F ≤ (K−1)2 and rank condition rank(Σ) = F. Global injectivity is proved rigorously by a three-step composition argument: (1) the linear parameterization θ ↦ Cθ is injective modulo the gauge; (2) the Sinkhorn map C ↦ π(C, a, b) is injective modulo the gauge (strict convexity of the dual → unique potentials → unique plan); (3) the normalization π ↦ Q is linear injective (ai > 0).

T2 (Sparsistency): The l1-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): Introduces the feature-moment map M(θ) = Φ⊤xθ. Local strong monotonicity is derived from the core spectral bound without compactness; global strong monotonicity 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 parameter µ ≥ π2min λmin(Σ)/ε2. Gradient descent converges monotonically to a local minimum; the empirical Hessian concentrates at the population Hessian at rate O(n−1/2).

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.

T1 (identifiability) is the foundation; T3 and T4 are both corollaries of the core spectral bound, 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 coupling (θ, ε) ↦ (cθ, cε).

The spectral proxy δxSSP = −(1/ε)PT Dπ PT δc reads as project to the tangent space → multiply elementwise by π → project back to the tangent space. The exact linearization is δx = −BHT−1B⊤δc, with Schur-complement form δx = −(1/ε)[Dπ − DπA⊤(ADπA⊤)−1ADπ]δc.

The score bound is sk(i, j) ≤ (πmax/(ε πmin))∥PT φk∥2, and the exponential lower bound for πmin is πmin(ε) ≥ amin bmin e−2∆C/ε, which cannot be replaced by a polynomial bound cεα.

The paper reports numerical verification across ten experiments (E1–E10): support recovery probability rising from 0.08 at n=200 to 0.96 at n=20000; condition-number divergence as ε→0 (rising by factor 3.8×104 from ε=0.5 to ε=0.02); the ε′-bias V-shape with empirical minimum at the generating ε=0.2; multiscale initialization raising success rate from 0.25 to 0.95; and misspecification projection residual 4.50 (vs. 10.63 for random models). The empirical Hölder exponent αeff ∈ 0.353, 0.266, 0.369, 0.298 across settings, and πmin scaling exponent αeff ≈ 19.9.

Improvement 1: Provably Identifiable Cost Learning in Inverse Problems

• What it does: An AI system can now learn the underlying cost function from observed optimal transport plans (e.g., matching data, routing logs) with a rigorous guarantee of uniqueness—up to a known gauge symmetry—provided the feature dimension satisfies F ≤ (K−1)2 and the feature covariance is full rank.

• Specific capability: Given empirical transition operators from entropic OT, the system outputs a cost parameter θ that is globally identifiable, eliminating the common failure mode of multiple equally good explanations. It can also detect when the data is insufficient (rank deficiency) and refuse to produce a spurious estimate.

Improvement 2: Sparse Support Recovery with Statistical Guarantees

• What it does: Adds an l1-penalized estimator that recovers the exact set of nonzero cost features (e.g., which pairwise attributes matter) with high probability, under explicit conditions: score concentration, irrepresentability, and a global selection condition.

• Specific capability: In high-dimensional settings (e.g., thousands of candidate features), the system selects only the truly relevant cost drivers (e.g., travel time, congestion, but not road color) with a sample complexity bound of O(∆2max log F / n). It can also report a confidence level for each selected feature, enabling safe automated feature discovery in logistics or economics.

Improvement 3: Well-Posed Inverse Mapping with Lipschitz Stability

• What it does: Guarantees that small perturbations in observed transport plans lead to bounded, predictable changes in the estimated cost parameter, with an explicit Lipschitz constant L ≤ ε∥Φ⊤Sa∥op / (πmin λmin(Σ)).

• Specific capability: An AI system can now quantify worst-case estimation error from noisy or partially observed OT data. For example, if the observation noise is 1%, the system can certify that the cost parameter error is at most 5% (given known constants), enabling deployment in safety-critical applications like medical image registration or autonomous vehicle trajectory matching.

Improvement 4: Monotone Convergent Optimization with Rate Guarantees

• What it does: Provides a gradient-descent algorithm for the cross-entropy objective that is guaranteed to converge monotonically to a local minimum, with a strong convexity parameter µ ≥ π2min λmin(Σ)/ε2.

• Specific capability: The system can train cost models from OT data without fear of divergence or slow oscillation. It can predict the number of iterations needed to reach a given accuracy (e.g., O(1/µ log(1/δ))), and it can automatically adapt the step size based on the spectral bound, making the training process reproducible and auditable.

Improvement 5: Robustness to Model Misspecification with Projection Consistency

• What it does: Under a compact parameter domain and a unique pseudo-true point, the estimator converges to the projection of the true (possibly non-OT) data-generating process onto the OT model class.

• Specific capability: When the real-world data does not exactly follow entropic OT (e.g., due to unmodeled noise or behavioral heterogeneity), the system still returns a meaningful answer: the closest OT model, with a quantified residual. It can also detect misspecification severity by comparing the residual to a threshold, alerting the user when the OT assumption is too strong.

Improvement 6: Spectral Proxy for Fast Sensitivity Analysis

• What it does: Replaces the expensive exact Hessian inverse in the Sinkhorn linearization with a spectrally equivalent proxy δxSSP = −(1/ε)PT Dπ PT δc, which preserves all singular-value bounds but is computationally trivial.

• Specific capability: An AI system can perform real-time what-if analysis on OT plans (e.g., how a 10% change in shipping cost affects the optimal matching) without recomputing the full Sinkhorn algorithm. This enables interactive decision-support tools for supply chain managers or urban planners, with a guaranteed error bound relative to the exact linearization.

Improvement 7: Data-Efficient Support Recovery via Multiscale Initialization

• What it does: Uses the paper’s multiscale ε-scheduling (start with large ε, refine to target ε) to raise support-recovery success from 0.25 to 0.95 in experiments.

• Specific capability: The system can automatically select a sequence of entropic regularization parameters, using the ε′-bias V-shape (empirically minimized at the true ε) to warm-start optimization. This reduces sample complexity by an order of magnitude (e.g., achieving 0.96 recovery at n=20,000 instead of requiring n>100,000), making the method practical for small datasets.

Improvement 8: Certified Condition-Number Monitoring

• What it does: Provides a closed-form bound on the Jacobian’s smallest singular value, σmin(Jθ) ≥ (πmin/(amax ε))√λmin(Σ), which diverges as ε→0 (verified empirically: factor 3.8×104 from ε=0.5 to 0.02).

• Specific capability: The system can predict when numerical instability will occur as ε decreases, and automatically stop or switch to a more stable parameterization (e.g., using the scale coupling (θ, ε) ↦ (cθ, cε) to renormalize). This prevents silent failures in production systems that must handle near-deterministic OT (ε→0).

We develop the statistical and algorithmic theory of inverse optimal transport (IOT) under the feature-parameterized cost C theta(i,j) = -theta T phi(i,j). The core technical contribution is the Sinkhorn linearization -- the implicit-function sensitivity of the entropic OT plan to the cost -- together with its spectral proxy, a formula that is spectrally exact yet geometrically transparent. The restricted Hessian on the tangent space satisfies the spectral sandwich (pi min/epsilon) I <= H T-1 <= (pi max/epsilon) I, yielding the single core bound sigma min >= (pi min/(a max epsilon)) sqrt(lambda min(Sigma)) that drives the entire theory. On this core we establish four theorems and one observation. T1 (identifiability): theta is globally injective on the quotient of the gauge kernel, with dimension bound F <= (K-1) squared. T2 (sparsistency): the l1-penalized estimator recovers the true support under irrepresentability and score concentration, with exponential failure probability. T3 (well-posedness): the feature-moment map M(theta) = Phi T x theta is strongly monotone, and the inverse is Lipschitz with constant L <= epsilon Phi T S a op / (pi min lambda min(Sigma)). T4 (convergence): local strong convexity with mu >= pi min squared lambda min(Sigma) / epsilon squared guarantees monotone gradient descent convergence. O5 (misspecification): the estimator converges to the OT-model projection of the truth; the Holder continuity of the projection map is assessed numerically, yielding setting-dependent empirical exponents alpha eff in (0,1).