Two histograms can have almost disjoint supports and still describe nearly the same object: imagine translating a narrow bump by one pixel. Optimal transport recognizes this—provided we already know which bins are close. But what if the bins are genes, words, or abstract features, with no obvious geometry?
With Geert-Jan Huizing and Laura Cantini, we studied a way to close the loop: compare samples using a distance between features, then compare features using the resulting distance between samples. A consistent pair is a pair of Wasserstein singular vectors [1]. The construction is unsupervised: unlike supervised ground-metric learning [2], it does not start from labels saying which samples ought to be similar.
One data matrix, two families of distributions
Let \(X\in\mathbb R_+^{n\times m}\) have nonzero rows and columns, with \(n,m\ge2\). Normalize each row to a probability vector \(a_i\in\Delta_m\), and each column to \(b_k\in\Delta_n\). Thus rows describe samples over features, while columns describe features over samples.
For a symmetric nonnegative cost matrix \(C\in\mathbb R^{m\times m}\) with zero diagonal, define the unrooted transport cost
\[ \operatorname{OT}_C(a,a')= \min_{P\ge0,\ P\mathbf1=a,\ P^\top\mathbf1=a'} \sum_{k,\ell}P_{k\ell}C_{k\ell}. \]Every feasible plan has total mass one. This is a finite linear program; no entropy or square root is included. For background on this formulation, see [4] and [5]. A feature cost produces a sample cost through
\[ \Phi_a(C)_{ij}=\operatorname{OT}_C(a_i,a_j). \]Interchange rows and columns to obtain \(\Phi_b(D)\in\mathbb R^{m\times m}\) from a sample cost \(D\). The self-consistency equations are
\[ \Phi_a(C)=\mu D,\qquad \Phi_b(D)=\lambda C, \qquad \lambda,\mu>0. \]One can fix scale by \(\|C\|_\infty=\|D\|_\infty=1\), where the norm is the entrywise maximum, not a matrix operator norm. These equations resemble left and right singular-vector equations, but the maps are nonlinear. The unknown vectors have become distance matrices.
Three simple properties do much of the work
Proposition 1 (order, scaling, and stability). For admissible costs and \(s\ge0\),
\[ \begin{aligned} C\le C'&\Longrightarrow\Phi_a(C)\le\Phi_a(C'),\\ \Phi_a(sC)&=s\Phi_a(C),\\ \|\Phi_a(C)-\Phi_a(C')\|_\infty&\le\|C-C'\|_\infty. \end{aligned} \]Proof: compare the objectives of the same transport plan
The first two claims follow directly by minimizing linear costs. For the third, if every entry changes by at most \(\delta\), the objective of every feasible plan changes by at most \(\delta\sum_{k,\ell}P_{k\ell}=\delta\). Compare an optimizer for either cost in the other problem. \(\square\)
If costs are strictly positive off the diagonal and the normalized rows are distinct, the output costs are also strictly positive off the diagonal. Order and homogeneity then give a projective stability estimate. Define
\[ d_H(C,C')= \max_{k\ne\ell}\log\frac{C_{k\ell}}{C'_{k\ell}} -\min_{k\ne\ell}\log\frac{C_{k\ell}}{C'_{k\ell}}. \]If \(\alpha C'\le C\le\beta C'\), the same sandwich passes through \(\Phi_a\). Taking the smallest possible ratio \(\beta/\alpha\) proves \(d_H(\Phi_a(C),\Phi_a(C'))\le d_H(C,C')\). In logarithmic coordinates, multiplicative homogeneity becomes additive homogeneity. This is the order geometry of nonlinear Perron–Frobenius theory [3].
Nonexpansiveness does not imply strict contraction, uniqueness, or convergence of every normalized iteration. Those conclusions need more information.
A regularizer that also makes existence transparent
The paper adds a simple term, different from entropic regularization:
\[ \Phi_a^\tau(C)_{ij} =\operatorname{OT}_C(a_i,a_j)+\tau\|C\|_\infty R(a_i-a_j), \qquad\tau>0, \]where \(R\) is a norm; define \(\Phi_b^\tau\) similarly. This term separates distinct histograms even when the ground cost has zero directions. Its sup-norm Lipschitz bound is \(1+\tau\max_{i,j}R(a_i-a_j)\), rather than the raw constant one.
Theorem 1 (regularized self-consistency). Assume the normalized rows are pairwise distinct and the normalized columns are pairwise distinct. For every \(\tau>0\), there are genuine metric matrices \(C,D\), positive off the diagonal, and positive numbers \(\lambda,\mu\), such that
\[ \Phi_a^\tau(C)=\mu D,\qquad \Phi_b^\tau(D)=\lambda C. \]Proof: normalize on a convex slice of the metric cone
First, transport preserves the triangle inequality. If \(C\) is a pseudometric, glue optimal plans from \(a\) to \(b\) and from \(b\) to \(d\). The pointwise inequality \(C_{k\ell}\le C_{kj}+C_{j\ell}\), integrated against their joint three-way law, gives
\[ \operatorname{OT}_C(a,d) \le\operatorname{OT}_C(a,b)+\operatorname{OT}_C(b,d). \]The term \(R(a-d)\) also obeys the triangle inequality. Thus both regularized transforms preserve pseudometrics, and their outputs are strictly positive on distinct histograms.
Next, choose the right normalization. Let \(\mathcal K_m\) be the set of symmetric pseudometric matrices with zero diagonal and \(\sum_{k<\ell}C_{k\ell}=1\). It is a nonempty compact convex set: all defining inequalities are linear. Use the analogous \(\mathcal K_n\) for \(D\).
Normalize each transformed matrix by the sum of its upper-triangular entries. The map
\[ (C,D)\longmapsto \left( \frac{\Phi_b^\tau(D)}{\sum_{k<\ell}\Phi_b^\tau(D)_{k\ell}}, \frac{\Phi_a^\tau(C)}{\sum_{i<j}\Phi_a^\tau(C)_{ij}} \right) \]is continuous and maps \(\mathcal K_m\times\mathcal K_n\) into itself. The denominators are positive by distinctness and \(\tau>0\). Brouwer’s theorem gives a fixed point. Undoing the normalization gives the two singular-vector equations, and the regularizer makes both metrics positive off the diagonal. Homogeneity then permits rescaling each to maximum one. \(\square\)
This proof uses a convex sum-normalized slice; the superficially more familiar constraint \(\|C\|_\infty=1\) is not itself a convex slice. The paper’s Theorem 2.3 gives regularized existence, and Proposition 2.6 proves uniqueness and linear convergence for sufficiently large regularization. We do not claim that every positive \(\tau\) is already in that contraction regime.
A power iteration, not a magic convergence theorem
The numerical update alternates the two transforms:
\[ C^{k+1}=\frac{\Phi_b^\tau(D^k)}{\|\Phi_b^\tau(D^k)\|_\infty}, \qquad D^{k+1}=\frac{\Phi_a^\tau(C^{k+1})}{\|\Phi_a^\tau(C^{k+1})\|_\infty}. \]Each cycle is a collection of transport problems. Starting inside the metric cone keeps every iterate there. At \(\tau=0\), the paper exhibits special cases with many self-consistent geometries; existence and uniqueness cannot be inferred from the raw nonexpansive estimate alone. At large \(\tau\), the histogram-norm term increasingly determines the learned geometry. These are useful checks on what the method can learn from the data.
Translations, scales, and an exact reference geometry
For periodic translated bumps, use \(X_{ik}=h((i-k)/N)\). With a translation-invariant initialization, the transforms stay circulant. A row of the learned cost matrix is then a cost profile versus circular displacement. The original paper explores this symmetry in its translated-histogram experiments [1]. A distance-like profile is plausible; a quadratic profile should be tested rather than assumed.
A second family varies both location and scale. Let \(h\) have mean zero and variance one, and set \(\rho_{m,s}(x)=s^{-1}h((x-m)/s)\), for \(s>0\). Its quantile is \(m+sQ_h\), so the one-dimensional quantile formula gives
\[ W_2^2(\rho_{m,s},\rho_{m',s'})=(m-m')^2+(s-s')^2. \]Proof: expand the quantile difference
Integrate \([(m-m')+(s-s')Q_h(u)]^2\) on \((0,1)\). The mixed term vanishes because \(\int Q_h=0\), and \(\int Q_h^2=1\). \(\square\)
The experiment compares learned sample costs to both this squared cost and its square root, after fitting one positive scale factor. A learned metric need not equal either. In particular, squared Euclidean cost generally violates the triangle inequality, while our iterates satisfy it. The paper’s debiased entropic variant has a separate large-regularization, PCA-related squared-cost limit; it is not the raw transform simulated here.
References
- Geert-Jan Huizing, Laura Cantini, and Gabriel Peyré. Unsupervised Ground Metric Learning Using Wasserstein Singular Vectors. ICML, PMLR 162, 9429–9443, 2022.
- Marco Cuturi and David Avis. Ground Metric Learning. Journal of Machine Learning Research 15, 533–564, 2014.
- Bas Lemmens and Roger Nussbaum. Nonlinear Perron–Frobenius Theory. Cambridge University Press, 2012.
- Gabriel Peyré and Marco Cuturi. Computational Optimal Transport. Foundations and Trends in Machine Learning 11(5–6), 355–607, 2019.
- Gabriel Peyré. Optimal Transport for Machine Learners. Book manuscript, 2025–2026.