Star-Separable Transport Maps in Semi-Discrete OT

• Star-separable transport maps are a property in semi-discrete optimal transport that partition the source domain into star-shaped Laguerre cells centered at discrete target points.
• The structure guarantees that each Laguerre cell contains its generator, facilitating accurate numerical integration and rapid convergence using Newton’s method.
• Leveraging this property yields significant computational improvements over traditional quasi-Newton and boundary discretization methods in multi-dimensional transport problems.

Star-separable transport maps are a structural property of solutions to semi-discrete optimal transport (OT) problems, characterized by a partitioning of the source domain into regions, each associated with a discrete target point. This property is established for cost functions that are positive combinations of $p$-norms, $1<p<\infty$, resulting in each region (also referred to as a Laguerre cell) being star-shaped with respect to its associated target point. The star-separable structure underlies efficient and accurate numerical algorithms for solving the Monge problem in two or more dimensions, with proven performance and convergence guarantees (Dieci et al., 2023).

1. Semi-Discrete Optimal Transport Formulation

The semi-discrete OT problem consists of finding a transport map $T:\Omega\to\{y_i\}_{i=1}^N$ between a continuous source measure and a discrete target measure that minimizes the transportation cost. The domain $\Omega\subset\mathbb{R}^d$ is compact, convex, with piecewise-smooth boundary and positive Lebesgue measure. The source measure $\mu$ is absolutely continuous: $d\mu(x)=\rho(x)dx,\;\rho\geq 0$. The target measure $\nu$ is discrete,

$$\nu = \sum_{i=1}^N m_i\,\delta_{y_i},\qquad m_i>0,\quad \sum_{i=1}^N m_i=1,$$

with $y_i\in \Omega^\circ$. The cost function takes the form

$$c(x,y) = \sum_{k=1}^m \alpha_k \|x-y\|_{p_k}, \qquad 1<p_k<\infty,\,\, \alpha_k>0,$$

which ensures positivity, symmetry, the triangle inequality, shift-invariance, and homogeneity. A partition of $\Omega$ into Laguerre cells $L_i(w)$ is generated via a weight vector $w\in\mathbb{R}^N$:

$$L_i(w) = \big\{x\in\Omega:\, c(x,y_i) - w_i \leq c(x,y_j) - w_j \; \forall j\neq i \big\}.$$

These Laguerre cells are disjoint and collectively exhaustive, i.e., $\Omega = \bigsqcup_{i=1}^N L_i(w)$.

2. Star-Shapedness of Laguerre Cells

The core result underpinning star-separable transport maps is the star-shapedness theorem: for any admissible cost $c(x,y)$ as above, each Laguerre cell $L_i(w)$ is star-shaped with respect to its generator $y_i$. Formally, for every $x\in L_i(w)$, the entire segment $[y_i, x]$ remains within $L_i(w)$. The argument is rooted in the metric properties of $c$, specifically the triangle inequality and shift-invariance.

Key lemmas justifying this structure include:

• Lemma A: each cell contains its generator, i.e., $y_i\in L_i(w)$ if $L_i(w)\neq\varnothing$.
• Lemma B: feasibility on boundaries demands $|w_i - w_j| \leq c(y_i, y_j)$ whenever $L_i(w)\cap L_j(w)\neq\emptyset$. The star-shapedness proof proceeds by contradiction using pathwise sampling along segments from $y_i$ to $x$ and the cost function’s convexity and triangle inequality, guaranteeing no exit from $L_i(w)$ along such rays.

3. Dual Potentials and Hessian Structure

The semi-discrete OT problem admits a dual in the weights $w$:

$$\Phi(w) = -\sum_{i=1}^N \int_{L_i(w)} [c(x,y_i) - w_i]\,d\mu(x) - \sum_{i=1}^N w_i m_i,$$

which is convex and continuously differentiable. The gradient is given by

$$\frac{\partial\Phi}{\partial w_i} = \mu(L_i(w)) - m_i,$$

so the optimal $w^*$ enforces that each cell has measure matching its target mass. The Hessian $H = \nabla^2 \Phi$ is, under mild regularity, symmetric and positive semi-definite with

$$H_{ij} = -\int_{L_i\cap L_j} \frac{\rho(x)}{\|\nabla_x c(x,y_i) - \nabla_x c(x,y_j)\|} d\sigma(x)\quad (i\neq j),\quad H_{ii} = -\sum_{j\neq i} H_{ij}.$$

Here $d\sigma$ is the $(d-1)$-dimensional surface element on the interface. The rank is $N-1$, and the nullspace is spanned by $e=(1,\dots,1)$.

4. Computational Methodology Leveraging Star-Separability

The star-shapedness yields a "star-separable" representation (Editor's term): each cell $L_i(w)$ can be parametrized by a single-valued radial function $r_i(\theta)$ in direction $\theta\in S^{d-1}$,

$$L_i(w) = \{ y_i + s\theta : \theta \in S^{d-1},\, 0 \leq s \leq r_i(\theta) \}.$$

This enables high-accuracy evaluation of integrals required for the gradient and Hessian:

• For the gradient:

$$\mu(L_i(w)) = \int_{S^{d-1}} \int_{0}^{r_i(\theta)} \rho(y_i + s\theta) s^{d-1} ds\, d\theta,$$

computed via adaptive 1-D quadrature (e.g., composite Simpson) over each directional slice.

• For the Hessian, off-diagonal entries reduce to surface integrals over facets $\partial L_i\cap\partial L_j$, split into sectors parameterized by $\theta$, evaluated again using 1-D adaptive quadrature. The star-shaped property guarantees $r_i(\theta)$ is single-valued, preventing ambiguity in this parametrization.

Newton’s method is applied to solve for $w$:

• The system $Hs=-\nabla\Phi$ is projected onto the subspace orthogonal to $e$; the Hessian’s singularity is inherent (total mass conservation).
• Step feasibility is checked: $|(w_i+s_i)-(w_j+s_j)|<c(y_i,y_j)$ for all $i\neq j$; steps are halved until feasibility.
• Descent in the merit function $g(w)=\| Q^\top \nabla\Phi \|^2$ can be enforced by line search for global convergence.

Initialization in 2D is feasible from $w=0$ (the Voronoi case), or by coarse-grid Newton–Raphson-like “auction” procedures.

5. Numerical Performance and Comparison

Extensive 2D numerical results are reported for the square $[0,1]^2$ using 2-norm and $p$-norm cost combinations. Empirical findings include:

• For $N=3$ to $N=10$ and both smooth/non-uniform source densities, Newton's method achieves convergence ($\|\nabla\Phi\|<10^{-8}$) in 2–6 steps.
• Radial (star-shaped) integration achieves area errors below machine tolerance ($10^{-12}$), significantly outperforming boundary discretization approaches.
• Analytic Hessian entries (surface integrals) are computed faster than finite-difference approximations of gradients.
• Comparisons with L-BFGS (quasi-Newton) and boundary subdivision methods reveal orders-of-magnitude improvements in speed and accuracy.
• Problems with $p\to 1$ or $p\to\infty$ remain tractable, though cell boundaries develop increased curvature, resulting in higher quadrature costs.

6. Scope and Extensions

The star-separable property for semi-discrete OT is valid for any cost function expressible as $\sum_k \alpha_k \|\cdot\|_{p_k}$ with $1 < p_k < \infty$. It provides both a conceptual simplification—reducing each cell to a radial structure—and a practical computational advantage, yielding rapid, globally convergent Newton solvers in low dimensions (notably 2D, with potential extension to 3D). The encapsulation of each Laguerre cell by a single-valued ray function is key for both geometric proofs and adaptive integration schemes, with direct impact on the efficiency and reliability of large-scale OT solvers (Dieci et al., 2023).