Theory Digest — Apr 30, 2026

1. Motivation: Understanding optimal representation learning requires principled theoretical foundations. The information bottleneck (IB) principle provides a framework for learning compressed representations that retain predictive information, but existing treatments lack geometric insight and practical non-variational algorithms.

2. Mathematical setup: Given random variables $(X,Y) \sim p(x,y)$ on measurable spaces $\mathcal{X} \times \mathcal{Y}$, learn encoder-decoder $X \to W \to \hat{Y}$ satisfying Markov chain $Y - X - W - \hat{Y}$. Define predictive manifold $\mathcal{M} = {p(Y

   x) : x \in \mathcal{X}}$ and probability simplex

   \[\mathcal{P}_K := \left\{p \in \mathbb{R}^K : p_i \geq 0, \sum_{i=1}^K p_i = 1\right\}\]

   where $K$ is effective predictive dimension. Assumptions:

   1. Standard regularity conditions for differential entropy

   2. Lipschitz regularity of $f^*: x \mapsto p(Y

      x)$ for continuous $Y$

   3. Finite effective dimension $K \leq d_X$ (intrinsic dimension of $\mathcal{X}$)

3. Toy example: Binary classification with $X$ uniform on ${1,2,3,4}$ and $Y \in {0,1}$ where $p(Y

   X=1) = p(Y

   X=2) = (0.1, 0.9)$ and $p(Y

   X=3) = p(Y

   X=4) = (0.9, 0.1)$. The predictive manifold $\mathcal{M} = {(0.1,0.9), (0.9,0.1)}$ contains two points in $\mathcal{P}_2$, and optimal encoder maps to these clusters.

4. Formal objective: Minimize encoding rate subject to predictive constraint:

   \[\min_{p(w|x)} I(X;W) \quad \text{s.t.} \quad I(X;Y|W) \leq \epsilon\]

The method develops a geometric framework connecting IB to rate-distortion theory with KL distortion:

\[d(x,w) := D_{KL}(p(y|x) \| p(y|w))\]

The optimal encoder has form:

\[p(w|x) = \frac{p(w)}{Z(x,\beta)} \exp(-\beta D_{KL}(p(y|x) \| p(y|w)))\]

Key steps:

1. Show optimal representation performs soft clustering of predictive manifold $\mathcal{M} \subset \mathcal{P}_K$
2. Derive transformation chain: flat Dirichlet $\to$ exponential $\to$ isotropic Gaussian
3. Connect to SIGReg via Cramér-Wold theorem for distributional regularization

4. Estimate encoder losses via minibatch marginals using CEB decomposition

   Applied to toy example: For binary case with two clusters, encoder assigns $x \in {1,2}$ to $w_1$ with $p(y

   w_1) = (0.1,0.9)$ and $x \in {3,4}$ to $w_2$ with $p(y

   w_2) = (0.9,0.1)$. As $\beta \to \infty$, assignment becomes deterministic. SIGReg enforces $W \sim N(0,I)$ in $\mathbb{R}^{2K}$ space with phase overhead $K \log 2\pi$ nats.

Step 1 — Prior work:

• “Information Bottleneck Method” (Tishby et al. 1999): introduced IB principle and self-consistency equations
• “β-VAE: Learning Basic Visual Concepts with a Constrained Variational Framework” (Higgins et al. 2017): variational approach to disentangled representations
• “Learning Representations by Maximizing Mutual Information” (Alemi et al. 2017): VIB with variational bounds

Step 2 — Delta: This paper provides geometric interpretation of IB as soft clustering of predictive manifold, derives exact Dirichlet-to-Gaussian transformation chain with quantified entropy overhead, and connects to SIGReg via non-variational estimation.

Step 3 — Theory-specific assessment:

• Main results are mostly predictable assemblies of known ingredients (rate-distortion equivalence, sufficient statistics, simplex geometry)
• Proof techniques are routine applications of information theory and differential entropy
• Bounds are not compared to lower bounds systematically
• The Dirichlet-to-Gaussian chain is elementary probability theory

Verdict: INCREMENTAL — solid synthesis of known results with clear geometric interpretation but no surprising theoretical advances.

Main proof strategy uses standard information-theoretic tools:

1. Rate-distortion equivalence via expectation manipulation:

   \[E[d(X,W)] = \sum_{x,w} p(x,w) \sum_y p(y|x) \log \frac{p(y|x)}{p(y|w)} = H(Y|W) - H(Y|X)\]

2. Optimal encoder via Lagrangian functional derivative:

   \[\frac{\delta}{\delta p(w|x)} L = \log p(w|x) + 1 + \beta \sum_y p(y|x) \log \frac{p(y|x)}{p(y|w)} - \lambda(x) = 0\]

3. Minimal sufficiency from zero-distortion constraint:

   \[I(X;Y|W^*) = 0 \Rightarrow d(x,w) = 0 \text{ a.s.} \Rightarrow p(y|w) = p(y|x)\]

4. Concavity of IB curve via standard convexity arguments for rate-distortion functions

5. Dirichlet-Gaussian chain using characteristic functions: If $X_1, \ldots, X_K \stackrel{iid}{\sim} \text{Exp}(1/2)$ and $X_k = Z_{k,1}^2 + Z_{k,2}^2$ with $Z_{k,j} \sim N(0,1)$, then normalized vector follows $\text{Dir}(1,\ldots,1)$

   The proofs are routine applications of classical information theory with no novel technical insights.

Experiments include:

1. Toy problems validating rate-distortion trade-offs on synthetic data
2. FashionMNIST classification comparing non-parametric estimator to standard variational approach

3. Verification of predicted IB curve shape and gap behavior

   Results show competitive performance of minibatch marginal estimation versus VIB, but experiments are limited in scope. The paper is primarily theoretical with empirical validation serving mainly as proof-of-concept rather than comprehensive evaluation.

Limitations:

1. Finite effective dimension assumption K ≤ d_X is TECHNICAL - needed for covering number bounds but likely removable with more sophisticated analysis
2. Lipschitz regularity requirement for continuous Y is RESTRICTIVE - significantly limits applicability to smooth predictive functions only
3. Minibatch estimation requires large batch sizes for accurate marginal estimation - NATURAL computational constraint

4. Phase entropy overhead K log 2π is non-diminishing - TECHNICAL limitation that affects rate accounting

   Open problems:

5. Characterize exact shape of IB curve beyond endpoints for general p(x,y)
6. Extend geometric framework to structured prediction tasks with combinatorial output spaces

Motivation: High-dimensional non-convex optimization problems like spike detection in random matrix models are ubiquitous in statistics and machine learning. Understanding when gradient-based methods succeed or fail is crucial, particularly whether initialization matters for convergence.

Mathematical setup: Consider the spiked Wigner matrix model:

\[M = G + \frac{\theta}{N} vv^T\]

where $G$ is an $N \times N$ GOE matrix with $G_{ij} = G_{ji} \sim \mathcal{N}(0, 1/N)$ for $i \neq j$ and $\mathcal{N}(0, 2/N)$ for $i = j$, $v \in S_N$ is uniformly random on the sphere $S_N = S^{N-1}(\sqrt{N})$, and $\theta > 0$ is the signal-to-noise ratio.

The Langevin dynamics on the sphere is:

\[dX_t = \nabla_{S_p} H(X_t) dt + \sqrt{\frac{2}{\beta}} dB_t\]

where

\[H(x) = \frac{1}{2}\langle x, Mx \rangle\]

and $\nabla_{S_p}$ is the spherical gradient. The Gibbs measure at inverse temperature $\beta > 0$ is:

\[\pi_\beta(A) \propto \int_A \exp(\beta H(σ)) dσ\]

Assumptions:

1. $G$ is a standard GOE matrix with independent Gaussian entries
2. $v$ is sampled uniformly from the sphere $S_N$ independently of $G$
3. $\beta = \alpha/\theta$ with $\alpha, \theta > 0$ fixed as $N \to \infty$

4. $\theta > \theta_0(\alpha)$ sufficiently large

   Toy example: When $N = 3$, $\theta = 2$, and $\alpha = 1.5$ (so $\beta = 0.75$), the matrix $M$ has approximate eigenvalues $\lambda_1 \approx 2.5$, $\lambda_2, \lambda_3 \approx 2$. The correlation $m_1(X_t) = \langle X_t, u_1 \rangle/\sqrt{N}$ evolves with drift away from zero when $

   m_1

   $ is small, creating metastable states at $\pm m_e$ where $m_e \approx \sqrt{1 - 2/3} \approx 0.58$.

   Formal objective: Determine the mixing time:

   \[t_{\text{mix}} = \inf\{t > 0 : \max_{x \in S_N} d_{TV}(P^t_x, \pi_\beta) \leq 1/4\}\]

The method analyzes mixing times by decomposing the problem based on initialization symmetry.

For symmetric mixing time (when initialization is invariant under $u_1 \mapsto -u_1$):

1. Show fast escape from equator: From any point with $

   m_1(X_0)

   = 0$, the correlation $m_1(X_t)$ hits threshold $m_3$ in time:

   \[T_{\text{hit}} = O\left(\frac{\log N}{(\lambda_1 - \lambda_2)(m_e^2 - m_3^2)}\right) = O(\log N)\]

2. Show retention of correlation: Once $

   m_1(X_t)

   > m_3$, the process stays above $m_2$ for sub-exponential time with high probability using comparison with mean-repelling OU process.

3. Apply local Bakry-Emery criterion: In regions ${

   m_1

   > m_{be}}$ where:

   \[m_{be} = \sqrt{\frac{\lambda_1 - \lambda_N - \frac{1}{\beta}}{2\lambda_1 - \lambda_2 - \lambda_N}}\]

   the condition $-\nabla^2_{S_p}H + \frac{1}{\beta}\text{Ric}_{S_p} \succeq \kappa_{be} I$ holds for $\kappa_{be} > 0$.

   For worst-case mixing time lower bound:

4. Compute free energy barrier: The probability of the equator ${m_1 = 0}$ under $\pi_\beta$ is:

   \[\pi_\beta(|m_1| \leq \epsilon) \leq \exp(-\Delta_{\alpha,\theta} N + C\epsilon N)\]

   where $\Delta_{\alpha,\theta}$ is the free energy difference between spiked and null models.

5. Use isoperimetric inequality: Lift to warped product manifold $S_N \times S^1$ and apply Cheeger’s inequality to bound spectral gap.

   Toy example application: For $\alpha = 1.5$, $\theta = 2$, we get $m_e \approx 0.58$, $m_{be} \approx 0.45$. Since $m_e > m_{be}$, fast symmetric mixing occurs in $O(\log 3) = O(1)$ time, while worst-case mixing takes $\exp(\Theta(3))$ time due to the $\pm u_1$ symmetry barrier.

Prior work:

1. Spiked matrix spectral analysis: Baik-Ben Arous-Péché (2005) established the BBP transition at $\theta_c = 1$ for eigenvalue/eigenvector recovery. Extensive follow-up work studied universality and fluctuations.

2. Langevin dynamics for spherical models: Belius-Kistler (2019), Ben Arous-Jagannath (2021) analyzed spherical spin glasses using dynamical mean field theory for $O(1)$ times, but required “warm starts” and didn’t address mixing times.

3. First-order methods for spiked models: Benaych-Georges et al. (2022) showed $O(1)$ recovery for $\theta = N^\delta$ (much stronger signal), but not information-theoretic thresholds.

   Delta: This paper provides the first rigorous analysis of Langevin mixing times at information-theoretic thresholds $\theta = O(1)$. Key innovations:

4. Initialization-dependent analysis: Distinguishes symmetric vs. worst-case mixing, showing the $\pm u_1$ symmetry is the sole bottleneck
5. Sharp metastability picture: Exact exponential rate $\Delta_{\alpha,\theta}$ for worst-case mixing via free energy computations

6. Technical innovation: Novel use of warped product manifolds to lift Langevin dynamics to apply isoperimetric inequalities

   Theory-specific assessment:

   • Surprising result: That symmetric initialization circumvents exponential bottleneck is non-obvious - one might expect other critical points to cause additional slowdowns
   • Proof technique: Genuinely new combination of: (i) SDE comparison with mean-repelling OU process, (ii) local Bakry-Emery conditions on caps, (iii) warped product lifting for isoperimetric bounds
   • Tightness: Lower bound from free energy barrier matches upper bound exactly, indicating optimality

   The threshold condition $\theta > \theta_0(\alpha)$ is technical but explicit. The warped product construction for applying Cheeger’s inequality to non-reversible dynamics is novel.

   Verdict: SIGNIFICANT — provides first rigorous mixing time analysis at information-theoretic thresholds with sharp constants, introducing new proof techniques for high-dimensional non-convex optimization.

The proof combines several sophisticated techniques:

1. SDE Comparison Theory: The key technical insight is comparing $m_1(X_t)^2$ to a mean-repelling Ornstein-Uhlenbeck process:

\[dY_t = Y_t dt + \frac{1}{\sqrt{N}} dW_t\]

Through the transformation $Y’_t = \sqrt{\frac{\beta C_E}{2}} \arcsin(m_t/C_E)$, the authors show:

\[dY_t'^2 \geq Y_t^2 dt + \text{noise terms}\]

when $

m_t

< m_3$, enabling hitting time bounds.

2. Local Bakry-Emery Analysis: On spherical caps ${

m_1

> m_{be}}$, the key inequality is:

\[-\nabla^2_{S_p}H + \frac{1}{\beta}\text{Ric}_{S_p} \succeq ((2\lambda_1 - \lambda_2 - \lambda_N)m^2 - (\lambda_1 - \lambda_N) - \frac{1}{\beta})I\]

Since $m_{be}^2 = \frac{\lambda_1 - \lambda_N - \frac{1}{\beta}}{2\lambda_1 - \lambda_2 - \lambda_N}$, this gives $\kappa_{be} > 0$ for $

m

> m_{be}$.

3. Free Energy Computations: Using cavity method techniques from spin glass theory, the equator probability satisfies:

\[\pi_\beta(|m_1| \leq \epsilon) \leq \exp(-\Delta_{\alpha,\theta} N + C\epsilon N)\]

where the free energy difference is:

\[\Delta_{\alpha,\theta} = \begin{cases}\] \[\frac{\alpha}{2} - \frac{(\alpha-1)^2}{4\theta^2} - \frac{1}{2}\log \alpha - \frac{1}{2} & 1 < \alpha \leq \theta \\\] \[\frac{\alpha}{2}(1-\frac{1}{\theta})^2 - \frac{1}{2}\log \theta - \frac{1}{4\theta^2} + \frac{1}{4} & \alpha > \theta\] \[\end{cases}\]

4. Warped Product Isoperimetric Analysis: The most novel technique lifts the Langevin dynamics to a warped product manifold $M = S_N \times S^1$ with metric:

\[g = g_{\text{round}}(x) + e^{2\beta H'(x)} d\theta d\theta\]

This enables the key identity: for $f \in C^\infty(S_N)$ and $\tilde{f} = f \circ p$,

\[\Delta_g \tilde{f} = Lf\]

where $L$ is the Langevin generator.

5. Metastable Analysis: The upper bound uses renewal theory - showing that from $\pi(·

m_1 \geq 0)$, hitting $m_1 = 0$ occurs with probability $e^{-\Delta_{\alpha,\theta}(1+o(1))N}$ per $O(\log N)$ equilibration period.

The technical innovation is handling the lifted geometry while controlling curvature bounds uniformly in $N$.

Purely theoretical.

Empirical validation would involve:

1. Numerical simulation of the Langevin SDE on spheres of varying dimension $N$
2. Measurement of mixing times from random vs. worst-case initializations
3. Verification of the phase transition at $\beta_c = 1/\theta$
4. Confirmation of the exponential rate $\Delta_{\alpha,\theta}$ for worst-case mixing

Limitations:

1. High signal-to-noise requirement: Conditions $\theta > \theta_{0,H}(\alpha)$ for fast mixing and $\theta > \theta_{0,L}(\alpha)$ for symmetric fast mixing are TECHNICAL - needed for proof technique but likely removable with more sophisticated analysis.

2. Gaussian disorder assumption: Results stated for GOE matrices, though authors claim extension to general Wigner matrices is straightforward - this is NATURAL.

3. Fixed $\alpha, \theta$ scaling: Analysis requires $\alpha, \theta = O(1)$ as $N \to \infty$. Behavior when these parameters scale with $N$ is not addressed - this is TECHNICAL.

4. Spherical prior: Assumes spike $v$ is uniformly random on sphere. Real applications may have different priors - this is NATURAL for the theoretical model.

5. Single spike: Only considers rank-1 spike. Multi-spike case has different dynamics - this is NATURAL as single-spike is the fundamental case.

   Open problems:

6. Optimal threshold: Determine the sharp threshold $\theta_0(\alpha)$ below which fast symmetric mixing fails, and whether it matches the information-theoretic threshold $\theta = 1$.

7. Scaling regimes: Analyze mixing times when $\theta, \alpha$ scale with $N$, particularly the critical window around $\theta = 1$ and the behavior for $\theta = 1 + N^{-\gamma}$.

1. Motivation: Optimal transport theory studies the minimum cost of transporting one probability measure to another. When using entropic regularization (the Sinkhorn divergence), the self-transportation cost becomes positive, creating bias. This paper studies when transport costs can be “debiased” by subtracting self-transportation terms.

2. Mathematical setup: Let $\mathcal{X}$ be a measurable space and $c: \mathcal{X} \times \mathcal{X} \to \mathbb{R} \cup {+\infty}$ be a symmetric cost function. Define the debiased cost:

   \[c_0(x,y) = c(x,y) - \frac{c(x,x)}{2} - \frac{c(y,y)}{2}\]

   Call $c$ debiasable if $c_0(x,y) \geq 0$ for all $x,y \in \mathcal{X}$.

   For $\epsilon \in [0,+\infty]$, the entropic optimal transport cost is:

   \[OT_c^\epsilon(\mu,\nu) = \inf_{\pi \in \Pi(\mu,\nu)} \int_{\mathcal{X} \times \mathcal{X}} c(x,y) d\pi(x,y) + \epsilon KL(\pi|\mu \otimes \nu)\]

   Assumptions (numbered list):

   1. $\mathcal{X}$ is a Polish space
   2. $c$ is lower semicontinuous and symmetric
   3. $\mu, \nu \in \mathcal{P}(\mathcal{X})$ (probability measures)

3. Toy example: When $\mathcal{X} = {1,2,3}$ and $c(1,3) = c(2,3) = 0$, $c(1,2) = 1$, all other costs zero. For $\mu = \frac{1}{2}\delta_1 + \frac{1}{2}\delta_2$ and $\nu = \delta_3$, we have $OT_c^\epsilon(\mu,\nu) = 0$ but $OT_c^\epsilon(\mu,\mu) > 0$, showing debiasability can fail.

4. Formal objective: Determine conditions under which $OT_c^\epsilon$ is debiasable:

   \[OT_c^\epsilon(\mu,\nu) - \frac{1}{2}OT_c^\epsilon(\mu,\mu) - \frac{1}{2}OT_c^\epsilon(\nu,\nu) \geq 0\]

The method establishes debiasability through inf-representations.

Key steps:

1. Characterize debiasable costs via inf-representation: any symmetric debiasable cost $c$ satisfies:

   \[c(x,y) = \inf_{z \in \mathcal{Z}} \psi(x,z) + \psi(y,z)\]

2. For classical OT ($\epsilon = 0$): Use probabilistic reformulation with random variables $T_1, T_2$:

   \[OT_c(\mu,\nu) = \inf_{(T_1)_\#\lambda = \mu, (T_2)_\#\lambda = \nu} \int c(T_1(w), T_2(w)) d\lambda(w)\]

3. For entropic OT ($\epsilon > 0$): - Case 1: Log-sum-exp costs. Define:

   \[c_{\epsilon,\lambda}(x,y) = -\epsilon \log \int_\mathcal{Z} e^{-\frac{\psi(x,z) + \psi(y,z)}{\epsilon}} d\lambda(z)\]

   - Case 2: Negative definite costs. Show $k(x,y) = e^{-c(x,y)/\epsilon}$ admits representation:
   

   \[k(x,y) = \int_\mathcal{Z} \rho_\epsilon(x,z)\rho_\epsilon(y,z) d\lambda(z)\]

4. For MMD limit ($\epsilon = +\infty$): Direct analysis using negative definiteness.

   Applied to toy example: For $c(x,y) =

   x-y

   ^2$ on $\mathbb{R}^d$, the inf-representation gives $\psi(x,z) = 2

   x-z

   ^2$, leading to:

   \[OT_c^\epsilon(\mu,\nu) = \inf_{\eta \in \mathcal{P}(\mathbb{R}^d)} 2W_2^2(\mu,\eta) + 2W_2^2(\nu,\eta) + \epsilon KL(\eta|\lambda)\]

Step 1 — Prior work:

• “Interpolating between Optimal Transport and MMD using Sinkhorn Divergences” (Feydy et al., 2019): First established debiasability for specific kernel-based costs on compact spaces
• “Sinkhorn divergences for unbalanced optimal transport” (Séjourné et al., 2019): Extended to unbalanced transport with similar assumptions

Step 2 — Delta: This paper provides:

• Unified characterization via inf-representations for all $\epsilon \in [0,+\infty]$
• Extension beyond compact spaces and universal kernels
• Novel decomposition formulas connecting entropic OT to Wasserstein barycenters
• Treatment of unbounded costs through bounded sets constraint

Step 3 — Theory-specific assessment:

• Main theorem is somewhat predictable as an extension of known compact/universal kernel results
• Proof techniques are largely routine, assembling inf-representation theory with entropic regularization
• The inf-representation characterization (Proposition 2.1) is not new, appearing in prior work by the authors
• Bounds are not compared to lower bounds (none appear to be known)

The connection to entropic barycenters and the $(ε, ε/2)$-scaling insight has some novelty, but the core debiasability results follow expected patterns.

Verdict: INCREMENTAL — solid extension of known results to broader settings, but techniques and main conclusions are predictable from prior compact/universal kernel theory.

Main proof strategy uses inf-representations and minimax techniques:

1. Classical OT (ε=0): Direct probabilistic argument via:

   \[\int c(T_1(w), T_2(w)) d\lambda(w) - \frac{1}{2}\int c(T_1(w), T_1(w)) d\lambda(w) - \frac{1}{2}\int c(T_2(w), T_2(w)) d\lambda(w) \geq 0\]

   following from debiasability of ground cost $c$.

2. Log-sum-exp costs: Key inequality from Cauchy-Schwarz:

   \[\left(\int_\mathcal{Z} e^{-\frac{\psi(x,z) + \psi(y,z)}{\epsilon}} d\lambda(z)\right)^2 \leq \int_\mathcal{Z} e^{-\frac{2\psi(x,z)}{\epsilon}} d\lambda(z) \int_\mathcal{Z} e^{-\frac{2\psi(y,z)}{\epsilon}} d\lambda(z)\]

3. Negative definite kernels: Utilize Gaussian series representation. For $c(x,y) =

   \phi(x) - \phi(y)

   ^2_{H_c} + s(x) + s(y)$:

   \[\rho_\epsilon(x,z) = \exp\left(\frac{2}{\epsilon}\sum_{i \in I} \phi_i(x)z_i - \frac{2}{\epsilon}\sum_{i \in I} \phi_i(x)^2 - \frac{s(x)}{\epsilon}\right)\]

4. Continuous kernels: Non-standard minimax argument exploiting weak compactness. Key technical insight is upper semicontinuity along maximizing sequences:

   \[\mathcal{L}(\mu', \nu', z) = \epsilon \int \log\left(\frac{d\mu'}{d\mu}\right) d\mu + \frac{1}{2}|k_{\mu'}|_k^2 + \frac{1}{2}|k_{\nu'}|_k^2 - (k_{\mu'}, k_{\nu'})_k\]

5. MMD limit: Direct computation using:

   \[-\frac{1}{2}\int c(x,y) d(\mu-\nu)(x) d(\mu-\nu)(y) = -\frac{1}{2}\int |\phi(x) - \phi(y)|^2_{H_c} d(\mu-\nu)(x) d(\mu-\nu)(y) \geq 0\]

Purely theoretical. The paper provides explicit computations for:

1. Squared Euclidean distance: Shows $OT_c^\epsilon(\delta_x, \delta_y) =

   x-y

   ^2 - \frac{d\epsilon}{2}\log(\pi\epsilon/4)$ with optimal auxiliary measure being Gaussian centered at midpoint $(x+y)/2$.

2. Finite discrete examples: Demonstrates failure of debiasability with 3-point space example.

   Empirical validation would require: (i) numerical verification of inf-representation formulas, (ii) comparison of computational efficiency between different representations, (iii) statistical properties of estimators based on these decompositions in practical applications.

Limitations:

1. TECHNICAL: Boundedness assumptions (e.g., $

   e^{c/\epsilon}

   _{L^1(\mu \otimes \nu)} \leq M$) needed for proof technique but likely removable with more sophisticated analysis

2. TECHNICAL: Restriction to Polish spaces for measurability - standard in optimal transport but could potentially be weakened

3. RESTRICTIVE: Negative definiteness requirement for general entropic case significantly narrows applicability beyond Euclidean/Gaussian settings

4. TECHNICAL: Integrally strictly positive definite condition for strict debiasability - needed for injectivity arguments but excludes some natural kernels

5. RESTRICTIVE: Compact space or bounded cost assumptions in Theorem 4.15 limit practical applicability to unbounded problems

   Open problems:

6. Computational complexity: What is the algorithmic complexity of computing the inf-representations compared to standard Sinkhorn iterations?

7. Statistical estimation: Develop finite-sample convergence rates for empirical versions of these debiased estimators, extending beyond the deterministic analysis presented.