Iterative Bregman Projections for Regularized Transportation Problems

Introduction

The theory of Optimal Transport (OT) defines a natural and useful geometry to compare measures supported on metric probability spaces. Its modern formulation as a linear program is due to Kantorovich . OT has recently found a flurry of applications in various fields such as computer vision , economy , computer graphics , image processing , astrophysics .

A major bottleneck that prevents the widespread of OT and its various generalizations is the lack of fast (possibly approximate) algorithms. Discrete optimal transport (i.e. computing transport between sums of Diracs) reduces to a finite dimensional linear program. When the mass of each Dirac is constant and the two measures have the same number $N$ of Dirac masses, this problem reduces to an optimal matching problem, for which dedicated discrete optimization methods exist , that roughly have $O(N^{3})$ complexity, which is still computationally too demanding for most applications. Another line of research, initiated by relies on dynamic formulations, which corresponds to computing the transport as a geodesic, and can be re-casted as a convex optimization problem. We refer to for an overview of several proximal optimization methods to tackle this problem. This requires adding an extra dimension (time variable along the geodesic) and is thus also computationally expensive. Semi-discrete optimal transport, i.e. optimal transport from a density to a weighted sum of dirac masses is a classical strategy for a generalised version of the Monge-Ampère equation, see for recent improvements of this approach. Finally and for the quadratic ground cost optimal transport, a direct Newton solver approach to the non-linear Monge Ampère equation can be used to compute density to density optimal transport. This holds under some regularity assumption on the densities and domain and hence on the transport map itself, see and for instance.

Entropic regularization.

A different approach consists in computing a regularized version of the OT problem. An interesting choice consists in penalizing the entropy of the joint coupling. This idea can be traced back to Schrodinger and can be related to the the so-called iterative proportional fitting procedure (IPFP) which has found numerous applications in the probability and statistics literature. We refer to for modern perspectives on this problem. Such a regularization also appears in the economy literature, where OT theory can be useful to predict flows of commodities or people in a market. In that context, regularizing the OT problem can also ensure the smoothness of such flows or facilitate inference in matching models . Such a regularization was also recently introduced in where it is shown that, in addition to favorable computational properties (parallelization, quadratic complexity) detailed below, such a regularization also yields a distance between histograms that can perform better in classification tasks than the usual OT distances.

The underlying idea of an entropic regularization is that entropy forces the solution to have a spread support, thus deviating from the fact that optimal couplings are sparse (i.e. supported on a graph of a transport plan solving Monge’s problem). A first impact of this regularization is that this non-sparsity of the solution helps to stabilize the computation. This can be related to the fact that entropic penalization defines a strongly convex program (as opposed to the initial OT problem) with a unique solution.

Another (even more important) advantage of this entropic regularized OT problem is that its solution is a diagonal scaling of $e^{-C}$, the element-wise exponential matrix of $-C$, where $C$ is the ground cost defining the transport (see Section 3.1 for more details). The solution to this diagonal scaling problem can be found efficiently through the IPFP iterative scheme . This algorithm was later studied in detail by Sinkhorn in and its convergence proof was extended to continuous measures in .

Entropic regularization of linear programs also shares some connection with interior point methods. These approaches make use of a $\log$-barrier function, which should be self-concordant to ensure a polynomial complexity for a given accuracy . Such a property does not hold for the entropic barrier, so it is not a competitive approach when it comes to approximating solutions of the original linear program by lowering the amount of regularization. As we advocate in this present paper, entropic regularization has however several other computational advantages (in particular because of its close connection with Kullback-Leibler projections), which makes it attractive when a slight amount of smoothing is acceptable in the computed approximation.

Optimization using the Kullback-Leibler Divergence.

When considering optimization over the simplex or the cone of positive vectors, it makes sense to replace the usual Euclidean metric by a divergence that quantifies with more relevance the difference between two vectors. Of particular interest for our work is the Kullback-Leibler (KL) divergence, since it is intimately related to entropic regularization. The simplest algorithmic block that can exploit such a divergence is the iterative projection on affine subsets of such cones under the KL divergence, which was introduced by Bregman . Computing the projection on the intersection of generic convex sets requires to replace iterative projections by more complicated algorithms, such as for instance Dykstra’s method . This algorithm is extended to Bregman divergences (such as KL) in and a proof of convergence is given in . For references in probability and statistics that also address the case of continuous distributions see . Note that several other proximal algorithms have been extended to this setting .

OT Barycenter.

The OT metric has been extended in many ways. A natural extension is to consider the barycenter between several distributions (the case of only 2 measures defining the usual transport). Such a barycenter is defined as the solution of a convex variational problem (a weighted sum of OT distances) over the space of measures, which is studied in details in . OT barycenters find applications for instance in statistics to define a mean empirical estimator from a family of observed histograms , or in machine learning to provide an extended definition of $k$-means clustering and compute average histograms-of-features under the OT metric.

Solving this variational problem is challenging. Two recent numerical works have addressed this problem: where a gradient descent on an entropic smoothing of OT distances is used and which is based on a dual formulation and tools from non-smooth optimization and computational geometry.

This barycenter problem can be extended to more complicated variational problem, such as for instance the Wasserstein propagation . Note that similar entropic regularization technics can be applied as well to this problem.

Multi-marginal transport.

The OT barycenter problem, as introduced in , is essentially equivalent to a multi-marginal optimal transport with quadratic cost as studied in . Multimarginal reformulations are not usually not tractable, since they involve an optimization problem whose size grows exponentially with respect to the number of marginals. Fortunately, the special structure of the OT barycenter problem leads to linear reformulations that are linear in the number of marginals, see and Section 3.2 below. There are however applications where the problem under study intrinsically has a multi-marginal structure, that cannot be factorized as barycenter computations.

Multi-marginal Optimal Transport, , is a natural extension of Optimal Transport with many potential fields of applications : Economics , Density Functionnal Theory in Quantum Chemistry . The first important instance of multi-marginal transport was probably Brenier’s generalised solutions of the Euler equations for incompressible fluids which are clearly described in his review paper . Note that entropic regularization of the multimarginal transport problem leads to a problem of multi-dimensional matrix scaling .

2 Contributions

In this paper, we present a unified framework to numerically solve entropic approximations of several generalized optimal transport problems. This framework corresponds to defining appropriate entropic penalizations of the initial linear programs. The key idea is then to interpret the corresponding problems as projections of some input Gibbs density on an intersection of convex sets according to the Kullback-Leibler divergence. This problem can then be solved efficiently using either Bregman iterative projection (for intersection of affine spaces) or a more general Dykstra-like algorithm—these being well-known first order non-smooth optimization schemes. We investigate in details the applications of these ideas to several generalized OT problems: barycenter (Section 3.2), tomographic reconstruction (Section 3.3), multi-marginal transport (Section 4), partial transport (Section 5.1) and capacity constrained transport (Section 5.2). The code implementing the methods presented in this article can be found onlinehttps://github.com/gpeyre/2014-SISC-BregmanOT.

3 Computational Speed

The goal of this paper is to present a new class of efficient methods to provide approximate solutions to linear program generalizing OT. It is however important to realize that these methods become numerically unstable when the regularization parameter (denoted $\varepsilon$ in the following) is small for two reasons: (i) since some of the quantities manipulated in the proposed algorithms have an order of $e^{-1/\varepsilon}$ (in particular Gibbs distributions denoted $\xi$ in the following) they become smaller than machine precision whenever the regularization $\varepsilon$ is small; (ii) more importantly, even if the first issue is taken care of by carrying out computations in the log domain, the convergence speed of iterative projection methods degrades significantly as $\varepsilon\rightarrow 0$. We observe therefore that these methods are competitive in a range where the regularization term $\varepsilon$ cannot be too small, and for which computed solutions exhibit a small amount of smoothing, see for instance Figure 1 for a visual illustration of this phenomenon. It is thus not the purpose of this article to compare these new methods with more traditional ones (such as interior points or simplex), because they do not target the same problem. Let us however single out the work of , that solves the regularized barycenter problem described in Section 3.2 using a gradient descent scheme. In all our numerical experiments, we found however that the iterative Bregman projection converge with substantially computational effort than this gradient descent and, because they rely on alternate projections, do not require adjusting gradient step-sizes and are thus easier to deploy.

4 Notations

The polytope of couplings between $(p,q)\in\Sigma_{N}^{2}$ is defined as

For a set $\mathcal{C}$, we denote $\iota_{\mathcal{C}}$ its indicator, that is

which is a concave function, where we used the convention $0\log(0)=0$.

With a slight abuse of notation, we extend these definitions for higher $d$-dimensional tensor arrays by replacing the sum over indices $(i,j)$ by sums over higher dimension indices.

Iterative Bregman Projections and Dykstra Algorithm

In this paper, we focus on regularized generalized OT problems that can be re-cast in the form

In the following, we extend the indexing of the sets by $L$-periodicity, so that they satisfy

One can then show that $\gamma^{(n)}$ converges towards the unique solution of (2),

2 Dykstra’s Algorithm

Dykstra’s algorithm starts by initializing

Recall here that $\odot$ and $\frac{\cdot}{\cdot}$ denotes entry-wise operations, see (1).

Dykstra algorithm converges to the solution of (2)

Entropic Regularization of Transport-like Problems

To illustrate the class of methods developed in this paper, we first review a classical approach to optimal transport approximation, that we recast in the language of Kullback-Leibler projections. This allows us to recover well known results, but in a framework that is easily generalizable.

Following many previous works (see Section 1.1 for details) we consider the following discrete regularized transport

The intuition underlying this regularization is that it enforces the optimal coupling $\gamma_{\varepsilon}^{\star}$ solution of (5) to be smoother as $\varepsilon$ increases. This regularization also yields favorable computational properties since problem (5) is $\varepsilon$-strongly convex. Its unique solution $\gamma_{\varepsilon}^{\star}$ can be obtained through elementary operations (matrix products, elementwise operations on matrices and vectors) as detailed below. If the optimal solution $\gamma^{\star}$ of the (original, non-regularized, i.e. $\varepsilon=0$) optimal transport problem is unique, then the optimal solution $\gamma_{\varepsilon}^{\star}$ of (5) converges to $\gamma^{\star}$ as $\varepsilon\rightarrow 0$. When other optimal solutions exist, $\gamma_{\varepsilon}^{\star}$ converges as $\varepsilon\rightarrow 0$ to that with the largest entropy among those, again denoted $\gamma^{\star}$. The convergence of minimizers of the regularized problem as $\varepsilon\to 0$ is actually exponential

where $\lambda$ and $M$ depend on $C$, $p$, $q$ and $N$ as shown by Cominetti and San Martin .

Problem (5) can be re-written as a projection

of $\xi$ according to the Kullback-Leibler divergence (here the exponential is computed component-wise).

The application of Bregman iterative projection (detailed in Section 2.1) to this splitting corresponds to the so-called IPFP/Sinkhorn algorithm (see Section 1.1 for bibliographical details).

The following well-known proposition details how to compute the relevant projections.

In plain words, the two projections in equation (7) normalize (with a multiplicative update) either the rows or columns of $\bar{\gamma}$ so that they have the desired row-marginal $p$ or column-marginal $q$.

An important feature of iterations (3), when combined with projections (7) is that the iterates $\gamma^{(n)}$ satisfy

This allows to implement this algorithm by only performing matrix-vector multiplications using a fixed matrix $\xi$, possibly in parallel if several OT are to be computed for several marginals sharing the same ground cost $C$ as shown in .

2 Optimal Transport Barycenters

We are given a set $(p_{k})_{k=1}^{K}$ of input marginals $p_{k}\in\Sigma_{N}$, and we wish to compute a weighted barycenter according to the Wasserstein metric. This problem finds many applications, as highlighted in Section 1.1.

Following , the general idea is to define the barycenter as a solution of a variational problem mimicking the definition of barycenters in Euclidean spaces. Given a set of normalized weights $\lambda\in\Sigma_{K}$, we consider the problem

It is easy to check that the Bregman iterative projection scheme can be applied to this setting by simply replacing $\operatorname{KL}$ by $\operatorname{KL}_{\lambda}$.

The $\operatorname{KL}_{\lambda}$ projection on $\mathcal{C}_{1}$ is computed as detailed in Proposition 1, since it is equal to the $\operatorname{KL}$ projection of each $\xi_{k}=\xi$ on a constraint of fixed marginal $p_{k}$. The $\operatorname{KL}_{\lambda}$ projection on $\mathcal{C}_{2}$ is computed as detailed in the following proposition.

where $\prod$ and $(\cdot)^{\lambda_{r}}$ should be understood as entry-wise operators.

Introducing the variable $p$ such that for all $k$, $\gamma_{k}\mathds{1}=p$, the first order conditions of the projection $P^{\tiny\operatorname{KL}_{\lambda}}_{\mathcal{C}_{2}}(\bar{\boldsymbol{\gamma}})$ states the existence of Lagrange multipliers $(u_{k})_{k}$ such that

Denoting $a_{k}=e^{-u_{k}}$, one has $\prod_{k}a_{k}=\mathds{1}$ and $\gamma_{k}=\operatorname{diag}(a_{k}^{1/\lambda_{k}})\bar{\gamma}_{k}$. Condition $\gamma_{k}\mathds{1}=p$ thus implies that

and condition $\prod_{k}a_{k}=\mathds{1}$ gives the desired value (10) for $p$. ∎

Note that when $K=2$, $(\lambda_{1},\lambda_{2})=(0,1)$, one retrieves exactly the IPFP/Sinkhorn algorithm to solve the entropic OT, as detailed in Section 3.1. Our novel scheme to compute barycenters should thus be understood as the natural generalization of this IPFP algorithm to barycenters.

Similarly as for Remark 1, one verifies that iterations (3) in the special case of problem (9) leads to iterates $\boldsymbol{\gamma}^{(n)}=(\gamma_{k}^{(n)})_{k}$ which satisfy, for each $k$

where $p^{(n)}$ is the current estimate of the barycenter, computed as

A nice feature of these iterations is that they can be computed in parallel for all $k$ using multiplications between the matrix $\xi$ and matrices storing $(u_{k}^{(n)})_{k}$ and $(v_{k}^{(n)})_{k}$ as columns.

Figure 2 shows an example of barycenters computation for $K=3$. The three vertices of the triangle show the input densities $(p_{1},p_{2},p_{3})$ which are uniform on binary shapes (diamond, annulus and square). The other points in the triangle display the results for the following values of $\lambda$

The computation is performed on an uniform 2-D grid of $N=256\times 256$ points in $^{2}$, and $\varepsilon=2/N$.

3 Partial Radon Inversion with OT Fidelity

The partial Radon transform (i.e. the computation of integrals of the data along parallel rays in a small limited set of directions) is a mathematical model for several scanning medical acquisition devices. This is an ill-posed linear operator, and inverting it while preventing noise and artifacts to blowup is of utmost importance for the targeted imaging applications. It is out of the scope of this paper to review the overwhelming literature on the topic of Radon inversion, and we refer to the book for an overview of classical approaches, and and the references therein for examples of state-of-the art methods.

The goal of this section is not to present a state of the art inversion scheme, but rather to show how the method recently introduced by can be solved using a simple iterative Bregman projection algorithm. We describe here the method in a fully discretized setting, where the Radon transform is implemented using a nearest neighbor interpolation.

We consider a square discretization grid of $N=N_{0}\times N_{0}$ pixels, indexed with

Given an angle $\theta$, we consider the following discrete lines in $\Omega_{N}$, $\forall\,(s_{1},s_{2})\in\Omega_{N}$,

where the mod $N_{0}$ is a modulo $N_{0}$ that maps the indices in the admissible range $\Omega_{N}$, hereby effectively implementing a convenient cyclic boundary condition.

The simplest way to perform a reconstructiong is to solve for the least squares estimate

This linear inverse does a poor job in the case of a small number $K$ of projections, since it exhibits reconstruction artifacts, as shown on Figure 3.

where $0<\lambda_{1}\leqslant 1$ is a weight that accounts for the degree of confidence in the template $g_{0}$, and $\lambda_{2}=1-\lambda_{1}$.

Here, $W_{2,\varepsilon}$ indicates the entropic Wasserstein distance (6) on the 2-D grid $\Omega_{N}$ where we defined the cost matrix $C=C^{2}$ as

and $W_{1,\varepsilon}$ indicate the entropic Wasserstein distance (6) on a 1-D periodic grid for the cost $C=C^{1}$

Similarly to the barycenter problem (9), we compute $f$ solving (13) as $f=\gamma\mathds{1}_{N}$ where $\gamma$ solves

(here the $\exp$ should be understood component-wise) and where we introduced

where the exponentiations are component-wise.

Figure 3 shows an example of application of the method for an image $f_{0}$ which is discretized on a grid of $N=80\times 80$ points, using $\varepsilon=2/N$ and $K=12$ Radon directions. There is no additional noise in the measurements, i.e. $w_{k}=0$ in (11). The template $g_{0}$ is a binary disk. These results show how using a large $\lambda_{1}$ recovers a result that is close to $g_{0}$, while using a smaller $\lambda$ introduces the geometric features of $f_{0}$ but also reconstruction artifacts. Note that the linear reconstruction $R^{+}(r)$ contains reconstruction artifacts.

Multi-marginal Optimal Transport

Multi-marginal optimal transport is a natural extension of optimal transport with many potential fields of applications, see Section 1.1.

The set of couplings between the marginals is

Similarly as (6), this problem can be re-cast as a KL projection

where the exponentiation is exponent-wise, and where

The Bregman projection on each of the convex $\mathcal{C}_{k}$ are again given by a simple normalisation as detailed in the following proposition.

For any $k$, denoting $\gamma=P^{\tiny\operatorname{KL}}_{\mathcal{C}_{k}}(\bar{\gamma})$, one has

It is thus possible to use the Bregman iterative projection detailed in Section 2.1 to compute the projection (16).

We now detail in the two following sections two typical cases of application of multi-marginal OT: grid-free barycenter computation (Section 4.2) and resolution of generalized Euler flow (Section 4.3).

2 Multi-marginal Barycenters

As shown in , the computation of barycenters of measures (thus the continuous analogous of (3.2)) can be computed by solving a multi-marginal transport problem.

An important point to note is that the measure barycenter (17) is in general composed of more than $N$ Diracs, and that these Diracs are not constrained to be on the discretization grid $(x_{i})_{i}$. In particular, the obtained result is different from the one obtained by solving (8), which computes a barycenter that lies on the same grid as the input measures. In some sense, formulation (17) is able to compute the “true” barycenter of measures, whereas (8) computes an approximation on a fixed grid, but the price to pay is the resolution of a high-dimensional multi-marginal program.

Figure 4 shows an histogram depiction of the measure $\mu_{\lambda}$ defined in (17), for the iso-barycenter (i.e. $\lambda_{k}=1/K$ for all $k$). It is computed by first solving (15) with the same three marginals used in Figure 2. The histogram $p\in\Sigma_{N}$ computed on a grid of $N=60\times 60$ points. Each $p_{i}$ is the total mass of $\mu_{\lambda}$ in the discretization square $S_{i}$ of size $1/\sqrt{N}\times 1/\sqrt{N}$, i.e. $p_{i}=\mu_{\lambda}(S_{i})$.

3 Generalized Euler Flows

Brenier proposed in a series of papers a relaxation of the Euler equation of incompressible fluids with constrained initial and final data. These data are conveniently expressed as a volume preserving map $\Xi$ of the domain. This relaxation can be understood as requiring the resolution of a multi-marginal transportation with an infinite number of marginals. Following (equation (21) section VII), when discretizing this problem with $K$ steps in time, one thus faces the resolution of a $K$ marginals OT problem.

We consider a fixed uniform discretization of $^{d}$ with points $(x_{i})_{i=1}^{N}$. The marginals are the uniform measure on this set (as discretization of the Lebesgue measure), i.e. $p_{k}=\mathds{1}/N$ for all $k$. The prescribed volume preserving maps $\Xi:^{d}\rightarrow^{d}$ is discretized using a permutation of the grid points, i.e. a discrete bijection $\sigma:\{1,\ldots,N\}\rightarrow\{1,\ldots,N\}$.

and the optimal coupling $\gamma$ solves (15).

It represents the evolution of a generalized flow of particles at time $t=\frac{k-1}{K-1}$. Note that in this setting, particles trajectories are non deterministic and their mass may split and spread across the domain.

Brenier’s numerical method is based on an approximation of the measure preserving map by a one to one permutation of the domain and the representation of the diffuse coupling therefore needs a large number of particles. Our resolution method is different and computes a space discretization of the coupling matrix and naturally encodes non-diffeomorphic volume preserving maps. The coupling $\gamma$ is an array of size $(N^{d})^{K}$ where the $d$-dimensional physical domain is discretized on $N^{d}$ points and we have $K$ times steps. However, as explained in Remark 4 below, because of the structure of the cost we only need to store and multiply $(N^{d})^{2}$ matrices.

and $D_{\alpha\beta}=|\!|x_{\alpha}-x_{\beta}|\!|^{2}$. We recall that all the marginals are equal to (a discretization of) the Lebesgue measure and $(x_{i})_{i}$ are discretized on the unit cube $^{d}$. As already noticed in Remark 1, the iterative Bregman projections (3) (always on the same Lebesgue marginal constraint) can be simplified as an IPFP iterative procedure. The optimal coupling $\gamma$ that solves (15) can be actually written as follows

and the IPFP procedure for $K$ marginals is

However, one can notice that the sum in (4) could be computationally onerous, but thanks to (19) we can rearrange it as

where $u^{k,(n)}$ is the $k^{\text{th}}$ vector at step $(n)$ and $B$ is the the product of $K$ smaller $N\times N$ matrices

Each iteration of the IPFP procedure therefore only involves $(2K)$ 2-coupling matrices multiplications and only requires storing $K$ vectors and the 2-coupling cost matrices $\xi^{0}$ and $\xi^{1}$. The computation of the 2-coupling maps (18) can be simplified with the same remark.

Figures 5, 6 and 7 show $T_{1,k}$ for three test cases in dimension $d=1$ proposed in . The computation is performed with a uniform discretization $(x_{i})_{i}$ of $$with$N=200$points,$\varepsilon=10^{-3}$and$K=16$. They agree with the solutions produced by Brenier and the mass spreading of the generalized flow is nicely captured by the 2 marginals couplings (18).

Transport Problems with Inequality Constraints

In this section, we consider transport problems with inequality constraints. Again we have to project for the KL divergence on the intersection of convex sets of nonnegative vectors.

The corresponding regularized problem reads

where the inequalities should be understood component-wise.

Similarly to (6), this is equivalent to computing the projection of $\xi=e^{-\frac{C}{\varepsilon}}$ on the intersection $\mathcal{C}_{1}\cap\mathcal{C}_{2}\cap\mathcal{C}_{3}$ of $K=3$ convex sets where

The following proposition shows that the KL projection onto those three sets can be obtained in closed form.

Since the considered sets $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ are convex but not affine, one thus needs to use Dykstra iterations (4) which are ensured to converge to the solution of (20).

If $\gamma^{\star}$ is the optimal solution of (20) and

are its marginals, then we define the active source $\mathcal{S}_{m}$ and the active target $\mathcal{T}_{m}$ regions as follow

where $\eta>0$ is a threshold we use to detect the region, namely the active region, where the transported mass is concentrated.

The continuous partial optimal transport problem has been studied in Caffarelli-McCann and Figalli . They show in particular that if there exists an hyperplane separating the support of the two marginals then the “active region” is separated from the “inactive region” by a free boundary which can be parameterized as a semi concave graph over the separating hyperplane. This can be observed on the test case presented in Figure 8. The computation is performed on an uniform 2D-grid of $N=256\times 256$ points in $^{2}$, $\varepsilon=10^{-3}$ and $m=0.7\min(\langle p,\,\mathds{1}\rangle,\langle q,\,\mathds{1}\rangle)$.

2 Capacity Constrained Transport

Korman and McCann proposed and studied in a variant of the classical OT problem when there is an upper bound on the coupling weights so as to capture transport capacity constraints.

where the inequalities should be understood component-wise.

This problem is equivalent to a KL projection problem of type (2) with $K=3$ convex sets and

The projection on $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ is given by Proposition 1. The projection on $\mathcal{C}_{3}$ is simply

where $R(x,y)=(x,-y)$ is the symmetry with respect to the second marginal axis and $W$ the optimal support of the saturated constraint $\mathcal{C}_{3}$.

3 Multi-Marginal Partial Transport

with $m\in[0,\min_{k}(\langle p_{k},\,\mathds{1}\rangle)]$

The KL projections on these convex sets are detailed in the following proposition.

For any $k=1,\ldots,K+1$, denoting $\gamma^{k}=P^{\tiny\operatorname{KL}}_{\mathcal{C}_{k}}(\bar{\gamma})$, one has

Once again the sets $\mathcal{C}_{k}$ are not affine, so one needs to use Dykstra iterations (4).

Figure 11 shows the results obtained when solving (25) with the same three marginals $(p_{1},p_{2},p_{3})$ used in Figure 2, using the cost

The computation is performed on an uniform 2D-grid of $N=60\times 60$ points in $^{2}$, $\varepsilon=0.005$ and $m=0.7\min_{k}(\langle p_{k},\,\mathds{1}\rangle)$.

Conclusion

In this paper, we have presented a unifying framework to approximate solutions of various OT-related linear programs through entropic regularization. This regularization enables the use of simple, yet powerful, iterative KL projection methods. While the entropy penalization is not a competitor with interior point methods when it comes to accurately solve the initial linear program, it produces fast approximations at the expense of an extra smoothing. It is thus a method of choice for many applications such as machine learning, image processing or economics.

Aknowledgements

We would like to thank Yann Brenier and Brendan Pass for stimulating discussions. The work of G. Peyré has been supported by the European Research Council (ERC project SIGMA-Vision). JD. Benamou, G. Carlier and L. Nenna gratefully acknowledge the support of the ANR, through the project ISOTACE (ANR-12-MONU-0013) and INRIA through the “action exploratoire” MOKAPLAN. M. Cuturi gratefully acknowledges the support of JSPS young researcher A grant 26700002 and the gift of a K40 card from the NVIDIA corporation.

References