Beyond entropic regularization: Debiased Gaussian estimators for discrete optimal transport and general linear programs

TL;DR

This work addresses the fundamental statistical challenge of inferring discrete OT plans with non-smooth dependence on marginals, which induces asymptotic bias for plug-in estimators and entropic-regularized methods. It introduces a dual-penalty regularization scheme in linear programs that yields an explicit, linear-in-$r$ bias term, enabling a Richardson-extrapolation-like debiasing via $\hat{\bm x}_n=2\bm x(r_n/2,\bm b_n)-\bm x(r_n,\bm b_n)$ that achieves Gaussian limits centered at the true solution. The authors develop a comprehensive theory for the trajectory and random perturbations of penalized LPs, establish bootstrap consistency for the debiased estimator, and demonstrate the approach through simulations and two real-data examples (Citi-Bike reallocation and colocalization analysis) showing improved inference over entropic regularization. The results extend beyond OT to standard-form LPs, providing a practical framework for asymptotically valid, data-driven inference in optimization-based problems. Overall, the paper delivers a principled, computationally tractable method for unbiased inference in discrete OT and broader LP settings, with broad applicability and strong empirical support.

Abstract

This work proposes new estimators for discrete optimal transport plans that enjoy Gaussian limits centered at the true solution. This behavior stands in stark contrast with the performance of existing estimators, including those based on entropic regularization, which are asymptotically biased and only satisfy a CLT centered at a regularized version of the population-level plan. We develop a new regularization approach based on a different class of penalty functions, which can be viewed as the duals of those previously considered in the literature. The key feature of these penalty schemes it that they give rise to preliminary estimates that are asymptotically linear in the penalization strength. Our final estimator is obtained by constructing an appropriate linear combination of two penalized solutions corresponding to two different tuning parameters so that the bias introduced by the penalization cancels out. Unlike classical debiasing procedures, therefore, our proposal entirely avoids the delicate problem of estimating and then subtracting the estimated bias term. Our proofs, which apply beyond the case of optimal transport, are based on a novel asymptotic analysis of penalization schemes for linear programs. As a corollary of our results, we obtain the consistency of the naive bootstrap for fully data-driven inference on the true optimal solution. Simulation results and two data analyses support strongly the benefits of our approach relative to existing techniques.

Figures (7)

• Figure 1: Examples of proper penalty functions (table on the left) and plot of these functions (figure on the right): log-barrier function (blue), inverse polynomial function with $\alpha=1$ (yellow), smoothed quadratic penalty function (green), exponential function (red).
• Figure 2: Comparison between the log and exponential penalty functions. Subfigures (a), (b) show the same experiment conducted with the logarithmic and exponential penalties, respectively, with 1000 independent replicates in each case. The first row shows plots of the mean-square error $\mathbb{E}\|\hat{\pi}_{n, r_0} - \pi^\star\|^2$ as the sample size $n$ and regularization parameter $r_0$ vary. The second shows the Kolmogorov--Smirnov distance between $\Delta w_{n, r_0}$ and the standard Gaussian in the same setting. In the last two rows, we fix $r_0=1$ and consider a small sample size $n=10^2$ and a large sample size $n=10^6$. The dashed line is a kernel density estimate of the density of $\Delta w_{n, r_0}$ and the solid line is the standard Gaussian density function as a reference. On the right of each density function plot is the corresponding Q--Q plot with a $45$-degree reference line in red.
• Figure 3: Large scale simulation: Experiments of optimal transportation between two independent Dirichlet distributions on equi-spaced $10\times10$ grid ($2000$ independent replicas in each case). With varies quantities of regularization parameter $r_0$ and sample size $n$, the left two plots shows the Mean Squared Error of the estimated plan ($\mathbb{E}\|\hat{\pi}_{n, r_0} - \pi^\star\|^2$), and the Kolmogorov–-Smirnov distance between $\Delta w_{n, r_0}$ and the standard Gaussian. The middle two plots show the density comparison between the standard Gaussian (solid lines) and the kernel density estimation of $\Delta w_{n, r_0}$ (dashed lines) with $r_0=1$ and $n=200$,$5000$ respectively, facilitated with the Q--Q plot on their right.
• Figure 4: Fast convergence with exponential penalty function in the $\bm{t}=\bm{s}$ case: The left two plots show the density comparison between the standard Gaussian (solid lines) and the kernel density estimation of $\sqrt{n}(\hat{w}_{n,1} - w^\star)$ (dashed lines) with $n=5000$ and $n=10^9$ respectively. The right log--log plot shows the rate of decay of $\mathbb{E}(\hat{w}_{n,1} - w^\star)^2$ with respect to $n$ in the large $n$ regime.
• Figure 5: Empirical bike transportation plan among all stations in Manhattan. The estimated quantities of bikes to be transported between station pairs are shown in purple lines. The two ends of a line are the locations of two stations. A thicker and darker line indicates that more bikes need to be transported between the two stations. Middle: overall plan. Left: Details of plan near Penn Station and Port Authority Bus Terminal. Upper right: Details of plan around Central Park. Lower right: Details of plan around Lower Manhattan and East Village.

Theorems & Definitions (55)

• Proposition 1.1
• Proposition 2.1
• Remark 1
• Lemma 3.1
• Definition 1: Decay rate of a function at $-\infty$
• Proposition 4.1
• Theorem 4.2: Convergence of trajectory of the penalized program
• Proposition 5.1
• Proposition 5.2
• Theorem 5.3