Flows, Scaling, and Entropy Revisited:
A Unified Perspective via Optimizing Joint Distributions

Consider a feasible flow for the OT problem in (1.1), i.e., a flow $f:E\to\mathbb{R}_{\geqslant 0}$ routing the supply distribution $\mu$ to the demand distribution $\nu$. This flow is naturally associated with a matrix $P\in\mathbb{R}_{\geqslant 0}^{n\times n}$ whose $ij$-th entry is the flow $f(i,j)$ on that edge. The netflow constraints on $f$ then simply amount to constraints on the row and column sums of $P$, namely $P1=\mu$ and $P^{T}1=\nu$. Thus, OT can be re-written as the LP

This decision set $\{P\in\operatorname*{\Delta_{n\times n}}:P1=\mu,\;P^{T}1=\nu\}$ is called the transportation polytope and can be equivalently viewed as the space of “couplings”—a.k.a., joint distributions $P(x,y)$ with first marginal $\mu$ and second marginal $\nu$.

MMC admits a similar formulation by taking an LP relaxation. Briefly, the idea is to re-write the objective in terms of matrices, as above, and then take the convex hull of the discrete decision set. Specifically, re-write the objective $\frac{1}{|\sigma|}\sum_{e\in\sigma}c(e)$ as $\langle P_{\sigma},C\rangle$ by associating to a cycle $\sigma$ the $n\times n$ matrix $P_{\sigma}$ with $ij$-th entry $1/|\sigma|$ if the edge $(i,j)\in\sigma$, and zero otherwise. By the Circulation Decomposition Theorem (see, e.g., [16, Problem 7.14]), the convex hull of the set $\{P_{\sigma}:\sigma\text{ cycle}\}$ of normalized cycles is the set $\{P\in\operatorname*{\Delta_{n\times n}}:P1=P^{T}1\}$ of normalized circulations, and so the LP relaxation of MMC is

and moreover this LP relaxation is exact. For details see, e.g., [3, Problem 5.47]. This decision set $\{P\in\operatorname*{\Delta_{n\times n}}:P1=P^{T}1\}$ can be equivalently viewed as the space of “self-couplings”—a.k.a., joint distributions $P(x,y)$ with symmetric marginals.

The most direct interpretation of the Sinkhorn algorithm is that it alternately projects³³3This projection is to the closest point in KL divergence rather than Euclidean distance. In the language of information geometry, this is an I-projection. the current matrix $P=XKY$ onto either the subspace $\{P\in\operatorname*{\Delta_{n\times n}}:P1=\mu\}$ with correct row marginals, or the subspace $\{P\in\operatorname*{\Delta_{n\times n}}:P^{T}1=\nu\}$ with correct column marginals. Note that when it corrects one constraint, it potentially violates the other. Nevertheless, the algorithm converges to the unique solution at the subspaces’ intersection, see Figure 1 for a cartoon illustration. The Osborne algorithm is analogous, except that it alternately projects $P=XKX^{-1}$ onto $n$ subspaces: the subspaces $\{P\in\operatorname*{\Delta_{n\times n}}:(P1)_{i}=(P^{T}1)_{i}\}$ defined by equal $i$-th row and column sums, for all $i\in[n]$. Here there is a choice for the order of subspaces to project onto; see [11] for a detailed discussion of this.

The Sinkhorn algorithm also admits an appealing dual interpretation. In the dual, entropic OT has $2n$ variables—the Lagrange multipliers $x,y\in\mathbb{R}^{n}$ respectively corresponding to the row and column marginal constraints in the primal. See Table 3, bottom right. Via the transformation $X_{ii}=e^{\eta x_{i}}$ and $Y_{ii}=e^{\eta y_{i}}$, these $2n$ dual variables are in correspondence with the $2n$ diagonal entries of the scaling matrices $X$ and $Y$ that the Sinkhorn algorithm seeks to find. One can verify that the Sinkhorn algorithm’s row update (Line 3 of Algorithm 1) is an exact block coordinate descent step that maximizes the dual regularized objective over all $x\in\mathbb{R}^{n}$ given that $y$ is fixed. Vice versa for the column update. See Figure 2 for a cartoon illustration. The Osborne algorithm is analogous, except that now there are only $n$ dual variables $x\in\mathbb{R}^{n}$ since there are only $n$ marginal constraints in the primal, see Table 4, bottom right. When the Osborne algorithm updates $X_{ii}=e^{\eta x_{i}}$ to equate the $i$-th entry of the row and column marginals, this corresponds to an exact coordinate descent step on $x_{i}$.

Each optimization problem in this note has been studied for many decades by many communities, which as already mentioned, has led to the development of many approaches besides the Sinkhorn and Osborne algorithms. Which algorithm is best? Currently there is no consensus. We suspect that the true answer is nuanced because different algorithms are often effective for different types of problem instances. For example, specialized combinatorial solvers blow competitors out of the water for small-to-medium problem sizes by easily achieving high accuracy solutions [26], whereas Sinkhorn and Osborne are the default algorithms in most numerical software packages for larger problems where moderate accuracy is acceptable. Both parameter regimes are important, but typically for different application domains. For example, high precision may be relevant for safety-critical or scientific computing applications. Whereas the latter regime is typically relevant for modern data-science applications of OT, since there is no need to solve beyond the inherent modeling error (OT is usually just a proxy loss in machine learning applications) and discretization error ($\mu,\nu$ are often thought of as samples from underlying true distributions). This latter regime is also relevant for pre-conditioning matrices before eigenvalue computation, since Matrix Scaling/Balancing optimize objectives that are just proxies for fast convergence of downstream eigenvalue algorithms [41, 46].

Comparisons between algorithms are further muddled by discrepancies between theory and practice. Sometimes algorithms are more effective in practice than our best theoretical guarantees suggest—is this because current analysis techniques are lacking, or because the input is an easy benchmark, or because of practical considerations not captured by a runtime theorem (e.g., the Sinkhorn algorithm is “embarrasingly parallelizable” and interfaces well with modern GPU hardware [25]). On the flipside, some algorithms with incredibly fast theoretical runtimes are less practical due to large constant or polylogarithmic factors hidden in the $\tilde{O}$ runtime. At least for now, this includes the elegant line of work (e.g., [54, 22, 4, 55]) based on the Laplacian paradigm, which very recently culminated in the incredible theoretical breakthrough [20] that solves the more general problem of minimum cost flow in time that is almost-linear in the input sparsity and polylogarithmic in $1/\varepsilon$. I am excited to see to what extent theory and practice are bridged in the upcoming years.

All algorithms discussed so far work for generic inputs. This robustness comes at an unavoidable $\Omega(n^{2})$ cost in runtime/memory from just reading/storing the input. This precludes scaling beyond $n$ in the several tens of thousands, say, on a laptop. For larger problems, it is essential to exploit “stucture” in the input. What stucture? This is a challenging question in itself because it is tied to the applications and pipelines relevant to practice. For OT, typically $\mu,\nu$ are distributions over $\mathbb{R}^{d}$ and the cost $C$ is given by pairwise distances, raised to some power $p\in[1,\infty)$. This is the $p$-Wasserstein distance, which plays an analogous role to the $\ell_{p}$ distance [56]. Then the $n\times n$ matrix $C$ is implicit from the $2n$ points in $\mu,\nu$. In low dimension $d$, this at least enables reading the input in $O(n)$ time—can OT also be solved in $O(n)$ time? A beautiful line of work in the computational geometry community has worked towards this goal, a recent breakthrough being $n\cdot(\varepsilon^{-1}\log n)^{O(d)}$ runtimes for $(1\raisebox{0.86108pt}{$\scriptstyle\pm$}\varepsilon)$ multiplicatively approximating OT for $p=1$ [47, 1]. We refer to those papers for a detailed account of this literature and the elegant ideas therein. Low-dimensional structure can also be exploited by the Sinkhorn algorithm. The basic idea is that the $n^{2}$ runtime arises only through multiplying the $n\times n$ kernel matrix $K=\exp[-\eta C]$ by a vector—a well-studied task in scientific computing related to Fast Multipole Methods [15]—and this can be done efficiently if the distributions lie on low-dimensional grids [53], geometric domains arising in computer graphics [53], manifolds [5], or algebraic varieties [13]. Preliminary numerics suggest that in these structured geometric settings, Sinkhorn can scale to millions of data points $n$ while maintaining reasonable accuracy [5]. Many questions remain open in this vein, for example graceful performance degradation in the dimension $d$, instance-optimality, and average-case complexity for “real-world inputs”.

This note focuses on optimizing joint distributions $P(x,y)$ with $k=2$ constrained marginals, and it is natural to ask about $k\geqslant 2$. These problems are called Multimarginal OT in the case of fixed marginals and linear objectives, and arise in diverse applications in fluid dynamics [17], barycentric averaging [2], graphical models [34], distributionally robust optimization [21], and much more; see the monographs [45, 44, 40]. The setting of symmetric marginals arises in quantum chemistry via Density Functional Theory [24, 19]. A central challenge in all these problems is that in general, it is intractable even to store a $k$-variate probability distribution, let alone solve for the optimal one. Indeed, a $k$-variate joint distribution in which each variable takes $n$ values is in correspondence with a $k$-order tensor that has $n^{k}$ entries—an astronomical number even for tiny $n,k=20$, say. As such, a near-linear runtime in the size of $P$ (the goal in this note for $k=2$) is effectively useless for large $k$, and it is essential to go beyond this by exploiting structure via implicit representations. See part II of the author’s thesis [6] and the papers upon which it is based [9, 7, 8, 10] for a systematic investigation of what structure enables $\mathrm{poly}(n,k)$ time algorithms, and for pointers to the extensive surrounding literature.

The original papers [43, 41] studied Matrix Balancing in the setting: given $K\in\mathbb{C}^{n\times n}$, find diagonal $X$ such that the $i$-th row and column of $XKX^{-1}$ have equal $\ell_{p}$ norm, for all $i\in[n]$. Definition 4 is equivalent for any finite $p$ (which suffices for this note) and elucidates the connection to optimizing distributions. See [11, §1.4] for details. For the case $p=\infty$, similar algorithms have been developed and were recently shown to converge in polynomial time in the breakthrough [51]. It would be interesting to reconcile the case $p=\infty$ as the analysis techniques there seem different and the connection to optimizing distributions seems unclear.

In this note, entropically regularized OT was motivated as a means to an end for computing OT efficiently. However, it has emerged as an interesting quantity in its own right and is now used in lieu of OT in many applications due to its better statistical and computational properties [25, 30, 38], as well as its ability to interface with complex deep learning architectures [45]. Understanding these improved properties is an active area of research bridging optimization, statistics, and applications. We are not aware of an analogous study of entropic MMC and believe this may be interesting.

I am grateful to Jonathan Niles-Weed, Pablo Parrilo, and Philippe Rigollet for a great many stimulating conversations about these topics. This work was partially supported by NSF Graduate Research Fellowship 1122374 and a TwoSigma PhD Fellowship.