Smooth and Sparse Optimal Transport

Introduction

Optimal transport (OT) distances (a.k.a. Wasserstein or earth mover’s distances) are a powerful computational tool to compare probability distributions and have recently found widespread use in machine learning (Cuturi, 2013; Solomon et al., 2014; Kusner et al., 2015; Courty et al., 2016; Arjovsky et al., 2017). While OT distances exhibit a unique ability to capture the geometry of the data, their application to large-scale problems has been largely hampered by their high computational cost. Indeed, computing OT distances involves a linear program, which takes super-cubic time in the data size to solve using state-of-the-art network-flow algorithms. Related to the Schrödinger problem (Schrödinger, 1931; Léonard, 2012), entropy-regularized OT distances have recently gained popularity due to their desirable properties (Cuturi, 2013). Their computation involves a comparatively easier differentiable and unconstrained convex optimization problem, which can be solved using the Sinkhorn algorithm (Sinkhorn and Knopp, 1967). Unlike unregularized OT distances, entropy-regularized OT distances are also differentiable w.r.t. their inputs, enabling their use as a loss function in a machine learning pipeline (Frogner et al., 2015; Rolet et al., 2016).

Despite its considerable merits, however, entropy-regularized OT has some limitations, such as introducing blurring in the optimal transportation plan. While this nuisance can be reduced by using small regularization, this requires a carefully engineered implementation, since the naive Sinkhorn algorithm is numerically unstable in this regime (Schmitzer, 2016). More importantly, the entropy term keeps the transportation plan strictly positive and therefore completely dense, unlike unregularized OT. This lack of sparsity can be problematic when the optimal transportation plan itself is of interest, e.g., in color transfer (Pitié et al., 2007), domain adaptation (Courty et al., 2016) and ecological inference (Muzellec et al., 2017). Sparsity in these applications is motivated by the principle of parsimony (simple solutions should be preferred) and by the enhanced interpretability of transportation plans.

Our contributions. This background motivates us to study regularization schemes that lead to smooth optimization problems (i.e., differentiable everywhere and with Lipschitz continuous gradient) while retaining the desirable property of sparse transportation plans. To do so, we make the following contributions.

We regularize the primal with an arbitrary strongly convex term and derive the corresponding smoothed dual and semi-dual. Our derivations abstract away regularization-specific terms in an intuitive way (§3). We show how incorporating squared $2$-norm and group-lasso regularizations within that framework leads to sparse solutions. This is illustrated in Figure 1 for squared $2$-norm regularization.

Next, we explore the opposite direction: replacing one or both of the primal marginal constraints with approximate smooth constraints. When using the squared Euclidean distance to approximate the constraints, we show that this can be interpreted as adding squared $2$-norm regularization to the dual (§4). As illustrated in Figure 1, that approach also produces sparse transportation plans.

For both directions, we bound the approximation error caused by regularizing the original OT problem. For the regularized primal, we show that the approximation error of squared $2$-norm regularization can be smaller than that of entropic regularization (§5). Finally, we showcase the proposed approaches empirically on the task of color transfer (§6).

An open-source Python implementation is available at https://github.com/mblondel/smooth-ot.

Background

If $f$ is strictly convex, then the supremum in (1) is uniquely achieved. Then, from Danskin’s theorem (1966), it is equal to the gradient of $f^{*}$:

The dual of a norm $\|\cdot\|$ is defined by $\|\bm{x}\|_{*}\coloneqq\sup_{\|\bm{y}\|\leq 1}~{}\bm{y}^{\top}\bm{x}.$ We say that a function is $\gamma$-smooth w.r.t. a norm $\|\cdot\|$ if it is differentiable everywhere and its gradient is $\gamma$-Lipschitz continuous w.r.t. that norm. Strong convexity plays a crucial role in this paper due to its well-known duality with smoothness: $f$ is $\gamma$-strongly convex w.r.t. a norm $\|\cdot\|$ if and only if $f^{*}$ is $\frac{1}{\gamma}$-smooth w.r.t. $\|\cdot\|_{*}$ (Kakade et al., 2012).

Optimal transport. We focus throughout this paper on OT between discrete probability distributions $\bm{a}\in\triangle^{m}$ and $\bm{b}\in\triangle^{n}$. Rather than performing a pointwise comparison of the distributions, OT distances compute the minimal effort, according to some ground cost, for moving the probability mass of one distribution to the other. The modern OT formulation, due to Kantorovich , is cast as a linear program (LP):

where $\mathcal{U}(\bm{a},\bm{b})$ is the transportation polytope

which is the so-called c-transform. Plugging it back into the dual, we get the “semi-dual”

For a recent and comprehensive survey of computational OT, see (Peyré and Cuturi, 2017).

Strong primal ↔↔\leftrightarrow Relaxed dual

We study in this section adding strongly convex regularization to the primal problem (3). We define

These assumptions are sufficient for (8) to be strongly convex w.r.t. $T\in\mathcal{U}(\bm{a},\bm{b})$. On first sight, solving (8) does not seem easier than (3). As we shall now see, the main benefit occurs when switching to the dual.

Let the (non-smooth) indicator function of the non-positive orthant be defined as

To define a smoothed version of $\delta$, we take the convex conjugate of $\Omega$, restricted to the non-negative orthant:

The optimal solution $T^{\star}$ of (8) can be recovered from $(\bm{\alpha}^{\star},\bm{\beta}^{\star})$ by $\bm{{{t}}}_{j}^{\star}=\nabla\delta_{\Omega}(\bm{\alpha}^{\star}+\beta_{j}^{\star}\bm{1}_{m}-\bm{{{c}}}_{j})\quad\forall j\in[n].$

For a proof, see Appendix A.1. Intuitively, the hard dual constraints $\alpha_{i}+\beta_{j}-{{c}}_{i,j}\leq 0~{}\forall i\in[m]~{}\forall j\in[n]$, which we can write $\sum_{j=1}^{n}\delta(\bm{\alpha}+\beta_{j}\bm{1}_{m}-\bm{{{c}}}_{j})$, are now relaxed with soft ones by substituting $\delta$ with $\delta_{\Omega}$.

2 Smoothed semi-dual formulation

We now derive the semi-dual of (8), i.e., the dual (11) with one of the two variables eliminated. Without loss of generality, we proceed to eliminate $\bm{\beta}$. To do so, we use the notion of smoothed max operator. Notice that

This is indeed true, since the supremum is always achieved at one of the simplex vertices. To define a smoothed max operator (Nesterov, 2005), we take the conjugate of $\Omega$, this time restricted to the simplex:

If $\Omega$ is $\gamma$-strongly convex over $\triangle^{m}\cap\operatorname*{dom}\Omega$, then $\text{max}_{\Omega}$ is $\frac{1}{\gamma}$-smooth and its gradient is defined by $\nabla\text{max}_{\Omega}(\bm{x})=\bm{y}^{\star}$, where $\bm{y}^{\star}$ is the supremum of (13). We next show that $\text{max}_{\Omega}$ plays a crucial role in expressing the conjugate of $\text{OT}_{\Omega}$.

Conjugate of $\text{OT}_{\Omega}$ w.r.t. its first argument

where $\Omega_{j}(\bm{y})\coloneqq\frac{1}{b_{j}}\Omega(b_{j}\bm{y})$.

A proof is given in Appendix A.2. With the conjugate, we can now easily express the semi-dual of (8), which involves a smooth optimization problem in $\bm{\alpha}$.

The optimal solution $T^{\star}$ of (8) can be recovered from $\bm{\alpha}^{\star}$ by $\bm{{{t}}}_{j}^{\star}=b_{j}\nabla\text{max}_{\Omega_{j}}(\bm{\alpha}^{\star}-\bm{{{c}}}_{j})\quad\forall j\in[n].$

Proof. $\text{OT}_{\Omega}(\bm{a},\bm{b})$ is a closed and convex function of $\bm{a}$. Therefore, $\text{OT}_{\Omega}(\bm{a},\bm{b})=\text{OT}_{\Omega}^{**}(\bm{a},\bm{b})$. $\square$

We can interpret this semi-dual as a variant of (7), where the max operator has been replaced with its smoothed counterpart, $\text{max}_{\Omega_{j}}$. Note that $\bm{\alpha}^{\star}$, as obtained by solving the smoothed dual (11) or semi-dual (15), is the gradient of $\text{OT}_{\Omega}(\bm{a},\bm{b})$ w.r.t. $\bm{a}$ when $\bm{\alpha}^{\star}$ is unique or a sub-gradient otherwise. This is useful when learning with $\text{OT}_{\Omega}$ as a loss, as done with entropic regularization in (Frogner et al., 2015).

Solving the optimization problems. The dual and semi-dual we derived are unconstrained, differentiable and concave optimization problems. They can therefore be solved using gradient-based algorithms, as long as we know how to compute $\nabla\delta_{\Omega}$ and $\nabla\text{max}_{\Omega}$. In our experiments, we use L-BFGS (Liu and Nocedal, 1989), for both the dual and semi-dual formulations.

3 Closed-form expressions

We derive in this section closed-form expressions for $\delta_{\Omega}$, $\text{max}_{\Omega}$ and their gradients for specific choices of $\Omega$.

Next, we present two choices of $\Omega$ that induce sparsity in transportation plans. The resulting dual and semi-dual expressions are new, to our knowledge.

Group lasso. Courty et al. (2016) recently proposed to use $\Omega(\bm{y})=\gamma(\sum_{i}y_{i}\log y_{i}+\mu\sum_{G\in\mathcal{G}}\|\bm{y}_{G}\|)$, where $\bm{y}_{G}$ denotes the subvector of $\bm{y}$ restricted to the set $G$, and showed that this regularization improves accuracy in the context of domain adaptation. Since $\Omega$ includes a negative entropy term, the same remarks as for negative entropy apply regarding the differentiability of $\delta_{\Omega}$ and smoothness of $\text{max}_{\Omega}$. Unfortunately, a closed-form solution is available for neither (10) nor (13). However, since the log keeps $\bm{y}$ in the strictly positive orthant and $\|\bm{y}_{G}\|$ is differentiable everywhere in that orthant, we can use any proximal gradient algorithm to solve these problems to arbitrary precision.

A drawback of this choice of $\Omega$, however, is that group sparsity is never truly achieved. To address this issue, we propose to use $\Omega(\bm{y})=\gamma(\frac{1}{2}\|\bm{y}\|^{2}+\mu\sum_{G\in\mathcal{G}}\|\bm{y}_{G}\|)$ instead. For that choice, $\delta_{\Omega}$ is smooth and is equal to

where $\bm{y}^{\star}$ decomposes over groups $G\in\mathcal{G}$ and equals

As noted in the context of group-sparse NMF (Kim et al., 2012), (17) admits a closed-form solution

where we defined $\bm{x}^{+}\coloneqq\frac{1}{\gamma}[\bm{x}]_{+}$. We have thus obtained an efficient way to compute exact gradients of $\delta_{\Omega}$, making it possible to solve the dual using gradient-based algorithms. In contrast, Courty et al. (2016) use a generalized conditional gradient algorithm whose iterations require expensive calls to Sinkhorn. Finally, because $\bm{{{t}}}_{j}^{\star}=\nabla\delta_{\Omega}(\bm{\alpha}^{\star}+\beta_{j}^{\star}\bm{1}_{m}-\bm{{{c}}}_{j})~{}\forall j\in[n]$, the obtained transportation plan will be truly group-sparse.

Relaxed primal ↔↔\leftrightarrow Strong dual

We now explore the opposite way to define smooth OT problems while retaining sparse transportation plans: replace marginal constraints in the primal with approximate constraints. When relaxing both marginal constraints, we define the next formulation:

where $\Phi(\bm{x},\bm{y})$ is a smooth divergence measure.

We may also relax only one of the marginal constraints:

where $\Phi(\bm{x},\bm{y})$ is defined as in Definition 2.

For both (19) and (20), the transportation plans will be typically sparse. As discussed in more details in §7, these formulations are similar to (Frogner et al., 2015; Chizat et al., 2016), with the key difference that we do not regularize $T$ with an entropic term. In addition, for $\Phi$, we propose to use $\Phi(\bm{x},\bm{y})=\frac{1}{2\gamma}\|\bm{x}-\bm{y}\|^{2}$, which is $\frac{1}{\gamma}$-smooth, while these works use a generalized Kullback-Leibler (KL) divergence, which is not smooth. Relaxing the marginal constraints is useful when normalizing input measures to unit mass is not suitable (Gramfort et al., 2015) or to allow for only partial displacement of mass. Relaxing only one of the two constraints is useful in color transfer (Rabin et al., 2014), where we would like all the probability mass of the source image to be accounted for but not necessarily for the reference image.

Dual interpretation. As we show in Appendix A.3, in the case $\Phi(\bm{x},\bm{y})=\frac{1}{2\gamma}\|\bm{x}-\bm{y}\|^{2}$, the dual of (19) can be interpreted as the original dual with additional squared $2$-norm regularization on the dual variables $\bm{\alpha}$ and $\bm{\beta}$. For the dual of (20), the additional regularization is on $\bm{\alpha}$ only (on the original dual or equivalently on the original semi-dual). For that choice of $\Phi$, the duals of (19) and (20) are strongly convex. The dual formulations are crucial to derive our bounds in §5.

Solving the optimization problems. While the relaxed and semi-relaxed primals (19) and (20) are still constrained problems, it is much easier to project on their constraint domain than on $\mathcal{U}(\bm{a},\bm{b})$. For the relaxed primal, in our experiments we use L-BFGS-B, a variant of L-BFGS suitable for box-constrained problems (Byrd et al., 1995). For the semi-relaxed primal, we use FISTA (Beck and Teboulle, 2009). Since the constraint domain of (20) has the structure of a Cartesian product $b_{1}\triangle^{m}\times\dots\times b_{n}\triangle^{m}$, we can easily project any $T$ on it by column-wise projection on the (scaled) simplex. Although not exlored in this paper, the block Frank-Wolfe algorithm (Lacoste-Julien et al., 2012) is also a good fit for the semi-relaxed primal.

Theoretical bounds

Convergence rates. The dual (11) is not smooth in $\bm{\alpha}$ and $\bm{\beta}$ when using entropic regularization but it is when using the squared $2$-norm, with constant upper-bounded by $\nicefrac{{n}}{{\gamma}}$ w.r.t. $\bm{\alpha}$ and $\nicefrac{{m}}{{\gamma}}$ w.r.t. $\bm{\beta}$. The semi-dual (15) is smooth for both regularizations, with the same constant of $\nicefrac{{1}}{{\gamma}}$, albeit not in the same norm. The relaxed and semi-relaxed primals (19) and (20) are both $\nicefrac{{1}}{{\gamma}}$-smooth when using $\Phi(\bm{x},\bm{y})=\frac{1}{2\gamma}\|\bm{x}-\bm{y}\|^{2}$. However, none of these problems are strongly convex. From standard convergence analysis of (projected) gradient descent for smooth but non-strongly convex problems, the number of iterations to reach an $\epsilon$-accurate solution w.r.t. the smoothed problems is $O(\nicefrac{{1}}{{\gamma\epsilon}})$ or $O(\nicefrac{{1}}{{\sqrt{\gamma\epsilon}}})$ with Nesterov acceleration.

Approximation error. Because the smoothed problems approach unregularized OT as $\gamma\to 0$, there is a trade-off between convergence rate w.r.t. the smoothed problem and approximation error w.r.t. unregularized OT. A question is then which smoothed formulations and which regularizations have better approximation error. Our first theorem bounds $\text{OT}_{\Omega}-\text{OT}$ in the case of entropic and squared $2$-norm regularization.

Approximation error of $\text{OT}_{\Omega}$

Let $\bm{a}\in\triangle^{m}$ and $\bm{b}\in\triangle^{n}$. Then,

Proof is given in Appendix A.4. Our result suggests that, for the same $\gamma$, the approximation error can often be smaller with squared $2$-norm than with entropic regularization. In particular, this is true whenever $\min\{H(\bm{a}),H(\bm{b})\}>\frac{1}{2}\min\left\{\|\bm{a}\|^{2},\|\bm{b}\|^{2}\right\}$, which is often the case in practice since $0\leq\min\{H(\bm{a}),H(\bm{b})\}\leq\min\{\log m,\log n\}$ while $0\leq\frac{1}{2}\min\left\{\|\bm{a}\|^{2},\|\bm{b}\|^{2}\right\}\leq\frac{1}{2}$. Our second theorem bounds $\text{OT}-\text{ROT}_{\Phi}$ and $\text{OT}-\widetilde{\text{ROT}}_{\Phi}$ when $\Phi$ is the squared Euclidean distance.

Approximation error of $\text{ROT}_{\Phi}$, $\widetilde{\text{ROT}}_{\Phi}$

Let $\bm{a}\in\triangle^{m}$, $\bm{b}\in\triangle^{n}$, $\Phi(\bm{x},\bm{y})=\frac{1}{2\gamma}\|\bm{x}-\bm{y}\|^{2}$. Then,

Proof is given in Appendix A.5. While the bound for $\widetilde{\text{ROT}}_{\Phi}$ is better than that of $\text{ROT}_{\Phi}$, both are worse than that of $\text{OT}_{\Omega}$, suggesting that the smoothed dual formulations are the way to go when low approximation error w.r.t. unregularized OT is important.

Experimental results

We showcase our formulations on color transfer, which is a classical OT application (Pitié et al., 2007). More experimental results are presented in Appendix C.

When $d(\bm{x},\bm{y})=\|\bm{x}-\bm{y}\|^{2}$, as used in our experiments, the above admits a closed-form solution: $\hat{\bm{x}}_{i}=\frac{\sum_{j=1}^{n}{{t}}_{i,j}\bm{y}_{j}}{\sum_{j=1}^{n}{{t}}_{i,j}}$. Finally, we use the new color $\hat{\bm{x}}_{i}$ for all pixels assigned to $\bm{x}_{i}$. The same process can be performed with respect to the $\bm{y}_{j}$, in order to transfer the colors in the other direction. We use two public domain images “fall foliage” by Bernard Spragg and “comunion” by Abel Maestro Garcia, and reduce the number of colors to $m=n=4096$. We compare smoothed dual approaches and (semi-)relaxed primal approaches. For the semi-relaxed primal, we also compared with $\Phi(\bm{x},\bm{y})=\frac{1}{\gamma}\text{KL}(\bm{x}||\bm{y})$, where $\text{KL}(\bm{x}||\bm{y})$ is the generalized KL divergence, $\bm{x}^{\top}\log\left(\frac{\bm{x}}{\bm{y}}\right)-\bm{x}^{\top}\bm{1}+\bm{y}^{\top}\bm{1}$. This choice is differentiable but not smooth. We ran the aforementioned solvers for up to $1000$ epochs.

Results. Our results are presented in Figure 2. All formulations clearly produced better results than unregularized OT. With the exception of the entropy-smoothed semi-dual formulation, all formulations produced extremely sparse transportation plans. The semi-relaxed primal formulation with $\Phi$ set to the squared Euclidean distance was the only one to produce colors with a darker tone.

2 Solver and objective comparison

We compared the smoothed dual and semi-dual when using squared $2$-norm regularization. In addition to L-BFGS on both objectives, we also compared with alternating minimization in the dual. As we show in Appendix B, exact block minimization w.r.t. $\bm{\alpha}$ and $\bm{\beta}$ can be carried out by projection onto the simplex.

Results. We ran the comparison using the same data as in §6.1. Results are indicated in Figure 3. When the problem is loosely regularized, we made two key findings: i) L-BFGS converges much faster in the semi-dual than in the dual, ii) alternating minimization converges extremely slowly. The reason for i) could be the better smoothness constant of the semi-dual (cf. §5). Since alternating minimization and the semi-dual have roughly the same cost per iteration (cf. Appendix B), the reason for ii) is not iteration cost but a convergence issue of alternating minimization. When using larger regularization, L-BFGS appears to converge slighly faster on the dual than on the semi-dual, which is likely thanks to its cheap-to-compute gradients.

3 Approximation error comparison

We compared empirically the approximation error of smoothed formulations w.r.t. unregularized OT according to four criteria: transportation plan error, marginal constraint error, value error and regularized value error (cf. Figure 4 for a precise definition). For the dual approaches, we solved the smoothed semi-dual objective (15), since, as we discussed in §5, it has the same smoothness constant of $\nicefrac{{1}}{{\gamma}}$ for both entropic and squared $2$-norm regularizations, implying similar convergence rates in theory. In addition, in the case of entropic regularization, the expressions of $\text{max}_{\Omega}$ and $\nabla\text{max}_{\Omega}$ are trivial to stabilize numerically using standard log-sum-exp implementation tricks.

Results. We ran the comparison using the same data as in §6.1. Results are indicated in Figure 4. For the transportation plan error and the (regularized) value error, entropic regularization required 100 times smaller $\gamma$ to achieve the same error. This confirms, as suggested by Theorem 1, that squared $2$-norm regularization is typically tighter. Unsurprisingly, the semi-relaxed primal was tighter than the relaxed primal in all four criteria. A runtime comparison of smoothed formulations is also important. However, a rigorous comparison would require carefully engineered implementations and is therefore left for future work.

Related work

Regularized OT. Problems similar to (8) for general $\Omega$ were considered in (Dessein et al., 2016). Their work focuses on strictly convex and differentiable $\Omega$ for which there exists an associated Bregman divergence. Following (Benamou et al., 2015), they show that (8) can then be reformulated as a Bregman projection onto the transportation polytope and solved using Dykstra’s algorithm . While Dykstra’s algorithm can be interpreted implicitly as a two-block alternating minimization scheme on the dual problem, neither the dual nor the semi-dual expressions were derived. These expressions allow us to make use of arbitrary solvers, including quasi-Newton ones like L-BFGS, which as we showed empirically, converge much faster on loosely regularized problems. Our framework can also accomodate non-differentiable regularizations for which there does not exist an associated Bregman divergence, such as those that include a group lasso term. Squared $2$-norm regularization was recently considered in (Li et al., 2016) as well as in (Essid and Solomon, 2017) but for a reformulation of the Wasserstein distance of order $1$ as a min cost flow problem along the edges of a graph.

Relaxed OT. There has been a large number of proposals to extend OT to unbalanced positive measures. Static formulations with approximate marginal constraints based on the KL divergence have been proposed in (Frogner et al., 2015; Chizat et al., 2016). The main difference with our work is that these formulations include an additional entropic regularization on $T$. While this entropic term enables a Sinkhorn-like algorithm, it also prevents from obtaining sparse $T$ and requires the tuning of an additional hyper-parameter. Relaxing only one of the two marginal constraints with an inequality was investigated for color transfer in (Rabin et al., 2014). Benamou (2003) considered an interpolation between OT and squared Euclidean distances:

While on first sight this looks quite different, this is in fact equivalent to our semi-relaxed primal formulation when $\Phi(\bm{x},\bm{y})=\frac{1}{2\gamma}\|\bm{x}-\bm{y}\|^{2}$ since (26) is equal to

However, the bounds in §5 are to our knowledge new. A similar formulation but with a group-lasso penalty on $T$ instead of $\frac{1}{2\gamma}\|T\bm{1}_{n}-\bm{a}\|^{2}$ was considered in the context of convex clustering (Carli et al., 2013).

Smoothed LPs. Smoothed linear programs have been investigated in other contexts. The two closest works to ours are (Meshi et al., 2015b) and (Meshi et al., 2015a), in which smoothed LP relaxations based on the squared $2$-norm are proposed for maximum a-posteriori inference. One innovation we make compared to these works is to abstract away the regularization by introducing the $\delta_{\Omega}$ and $\text{max}_{\Omega}$ functions.

Conclusion

We proposed in this paper to regularize both the primal and dual OT formulations with a strongly convex term, and showed that this corresponds to relaxing the dual and primal constraints with smooth approximations. There are several important avenues for future work. The conjugate expression (14) should be useful for barycenter computation (Cuturi and Peyré, 2016) or dictionary learning (Rolet et al., 2016) with squared $2$-norm instead of entropic regularization. On the theoretical side, while we provided convergence guarantees w.r.t. the OT distance value as the regularization vanishes, which suggested the advantage of squared $2$-norm regularization, it would also be important to study the convergence w.r.t. the transportation plan, as was done for entropic regularization by Cominetti and San Martín (1994). Finally, studying optimization algorithms that can cope with large-scale data is important. We believe SAGA (Defazio et al., 2014) is a good candidate since it is stochastic, supports proximity operators, is adaptive to non-strongly convex problems and can be parallelized (Leblond et al., 2017).

Acknowledgements

We thank Arthur Mensch and the anonymous reviewers for constructive comments.

References

Appendix A Proofs

We now add Lagrange multipliers for the two equality constraints but keep the constraint $T\geq 0$ explicitly:

Since (29) is a convex optimization problem with only linear equality and inequality constraints, Slater’s conditions reduce to feasibility (Boyd and Vandenberghe, 2004, §5.2.3) and hence strong duality holds:

Finally, plugging the expression of (10) gives the claimed result.

A.2 Derivation of the convex conjugate

The convex conjugate of $\text{OT}_{\Omega}(\bm{a},\bm{b})$ w.r.t. the first argument is

Following a similar argument as (Cuturi and Peyré, 2016, Theorem 2.4), we have

Notice that this is an easier optimization problem than (8), since there are equality constraints only in one direction. Cuturi and Peyré (2016) showed that this optimization problem admits a closed form in the case of entropic regularization. Here, we show how to compute $\text{OT}_{\Omega}^{*}$ for any strongly-convex regularization.

The problem clearly decomposes over columns and we can rewrite it as

where we defined $\Omega_{j}(\bm{y})\coloneqq\frac{1}{b_{j}}\Omega(b_{j}\bm{y})$ and where $\text{max}_{\Omega}$ is defined in (13).

A.3 Expression of the strongly-convex duals

Using a similar derivation as before, we obtain the duals of (19) and (20).

where $\Phi^{*}$ is the conjugate of $\Phi$ in the first argument.

The duals are strongly convex if $\Phi$ is smooth. When $\Phi(\bm{x},\bm{y})=\frac{1}{2\gamma}\|\bm{x}-\bm{y}\|^{2}$, $\Phi^{*}(-\bm{\alpha},\bm{a})=\frac{\gamma}{2}\|\bm{\alpha}\|^{2}-\bm{\alpha}^{\top}\bm{a}$. Plugging that expression in the above, we get

This corresponds to the original dual and semi-dual with squared $2$-norm regularization on the variables.

A.4 Proof of Theorem 1

Before proving the theorem, we introduce the next two lemmas, which bound the regularization value achieved by any transportation plan.

Bounding the entropy of a transportation plan

Let $H(\bm{a})\coloneqq-\sum_{i}a_{i}\log a_{i}$ and $H(T)\coloneqq-\sum_{i,j}{{t}}_{i,j}\log{{t}}_{i,j}$ be the joint entropy. Let $\bm{a}\in\triangle^{m}$, $\bm{b}\in\triangle^{n}$ and $T\in\mathcal{U}(\bm{a},\bm{b})$. Then,

Proof. See, for instance, (Cover and Thomas, 2006).

Together with $0\leq H(\bm{a})\leq\log m$ and $0\leq H(\bm{b})\leq\log n$, this provides lower and upper bounds for the entropy of a transportation plan. As noted in (Cuturi, 2013), the upper bound is tight since

Bounding the squared $2$-norm of a transportation plan

Let $\bm{a}\in\triangle^{m}$, $\bm{b}\in\triangle^{n}$ and $T\in\mathcal{U}(\bm{a},\bm{b})$. Then,

Proof. The tightest lower bound is given by $\displaystyle{\min_{T\in\mathcal{U}(\bm{a},\bm{b})}}\|T\|^{2}$. An exact iterative algorithm was proposed in (Calvillo and Romero, 2016) to solve this problem. However, since we are interested in an explicit formula, we consider instead the lower bound $\displaystyle{\min_{\begin{subarray}{c}T\bm{1}_{n}=\bm{a}\\ T^{\top}\bm{1}_{m}=\bm{b}\end{subarray}}}\|T\|^{2}$ (i.e., we ignore the non-negativity constraint). It is known (Romero, 1990) that the minimum is achieved at ${{t}}_{i,j}=\frac{{a}_{i}}{n}+\frac{{b}_{j}}{m}-\frac{1}{mn}$, hence our lower bound. For the upper bound, we have

We can do the same with $\bm{b}\in\triangle^{n}$ to obtain $\|T\|^{2}\leq\|\bm{b}\|^{2}$, yielding the claimed result. $\square$

Together with $0\leq\|\bm{a}\|^{2}\leq 1$ and $0\leq\|\bm{b}\|^{2}\leq 1$, this provides lower and upper bounds for the squared $2$-norm of a transportation plan.

Proof of the theorem. Let $T^{\star}$ and $T^{\star}_{\Omega}$ be optimal solutions of (3) and (8), respectively. Then,

Using $T^{\star},T^{\star}_{\Omega}\in\mathcal{U}(\bm{a},\bm{b})$ together with Lemma 46 and Lemma 48 gives the claimed results.

A.5 Proof of Theorem 24

To prove the theorem, we first need the following two lemmas.

Bounding the $1$-norm of $\bm{\alpha}$ and $\bm{\beta}$ for $(\bm{\alpha},\bm{\beta})\in\mathcal{P}(C)$

Let $\bm{\alpha},\bm{\beta}\in\mathcal{P}(C)$ with extra constraints $\bm{\alpha}^{\top}\bm{1}_{m}=0$ and $\bm{\alpha}^{\top}\bm{a}+\bm{\beta}^{\top}\bm{b}\geq 0$, where $\bm{a}\in\triangle^{m}$ and $\bm{b}\in\triangle^{n}$. Then,

Proof. The proof technique is inspired by (Meshi et al., 2012, Supplementary material Lemma 1.2).

Our goal is to upper bound the following objective

with a constant that does not depend on $\bm{r}$ and $\bm{s}$. We call the above the dual problem. Its Lagrangian is

By weak duality, any feasible primal point provides an upper bound of the dual problem. We start by choosing $\mu=\frac{1}{m}(\sum_{j}s_{j}-\sum_{i}r_{i})$ so that $\sum_{i,j}{{t}}_{i,j}$ provides the same values w.r.t. the last two constraints. Next, we choose

which ensures the non-negativity of $\nu\bm{a}+\bm{r}+\mu\mathbf{1}_{m}$ and $\nu\bm{b}+\bm{s}$ regardless of $\bm{r}$ and $\bm{s}$. It follows that the transportation plan $T$ defined by

is feasible. We finally bound the objective, $\langle T,C\rangle\leq\left|\!\left|C\right|\!\right|_{\infty}\sum_{i,j}{{t}}_{i,j}\leq\left|\!\left|C\right|\!\right|_{\infty}(\nu+n)$. $\square$

Bounding the $1$-norm of $\bm{\alpha}$ for $(\bm{\alpha},\cdot)\in\mathcal{P}(C)$

Let $\bm{\alpha},\bm{\beta}\in\mathcal{P}(C)$ with extra constraints $\sum_{i=1}^{m}\alpha_{i}=0$ and $\alpha^{\top}\bm{a}+\bm{\beta}^{\top}\bm{b}\geq 0$, where $\bm{a}\in\triangle^{m}$ and $\bm{b}\in\triangle^{n}$. Then,

Proof. Similarly as before, our goal is to upper bound

with a constant which does not depend on $\bm{r}$. The corresponding primal is

By weak duality, any feasible primal point gives us an upper bound. We start by choosing $\mu=\frac{1}{m}\sum_{i}r_{i}$ so that $\sum_{ij}{{t}}_{i,j}$ provides the same values w.r.t. the last two constraints. Next, we choose, $\nu=\underset{i}{\max}\frac{2}{a_{i}}$, which ensures the non-negativity of $\nu\bm{a}+\bm{r}+\mu\mathbf{1}_{m}$ ($\nu\bm{b}\geq 0$ is also satisfied since $\nu\geq 0$) which appears in the r.h.s. of the second constraint, independently of $\bm{r}$. It follows that the transportation plan $T$ defined by

is feasible. We finally bound the objective

Proof of the theorem. We begin by deriving the bound for the relaxed primal. Let $(\bm{\alpha}^{\star},\bm{\beta}^{\star})$ and $(\bm{\alpha}^{\star}_{\Phi},\bm{\beta}^{\star}_{\Phi})$ be optimal solutions of (5) and (43), respectively. Since $(\bm{\alpha}^{\star}_{\Phi})^{\top}\bm{a}+(\bm{\beta}^{\star}_{\Phi})^{\top}\bm{b}\leq(\bm{\alpha}^{\star})^{\top}\bm{a}+(\bm{\beta}^{\star})^{\top}\bm{b}$, we have

Taking the square of this bound and plugging the result in (72) gives the claimed result. Applying the same reasoning with Lemma 65 gives the claimed result for the semi-relaxed primal.

Appendix B Alternating minimization with exact block updates

General case. Let $\bm{\beta}(\bm{\alpha})$ be an optimal solution of (11) given $\bm{\alpha}$ fixed, and similarly for $\bm{\alpha}(\bm{\beta})$. From the first-order optimality conditions,

and similarly for $\bm{\alpha}$ given $\bm{\beta}$ fixed. Solving these equations is non-trivial in general. However, because

Entropic regularization. It is easy to verify that (75) is satisfied with

and similarly for $\bm{\alpha}(\bm{\beta})$. These updates recover the iterates of the Sinkhorn algorithm (Cuturi, 2013).

Squared $2$-norm regularization. Plugging the expression of $\nabla\delta_{\Omega}$ in (75), we get that $\bm{\beta}(\bm{\alpha})$ must satisfy

Close inspection shows that it is exactly the same optimality condition as the Euclidean projection onto the simplex $\displaystyle{\operatorname*{argmin}_{\bm{y}\in\triangle^{m}}}\|\bm{y}-\bm{x}\|^{2}$ must satisfy, with $\bm{x}=\frac{\bm{\alpha}-\bm{{{c}}}_{j}}{\gamma b_{j}}$. Let $x_{}\geq\dots\geq x_{[m]}$ be the values of $\bm{x}$ in sorted order. Following (Michelot, 1986; Duchi et al., 2008), if we let

then $\bm{y}^{\star}$ is exactly achieved at $[\bm{x}+\frac{\beta_{j}(\bm{\alpha})}{\gamma b_{j}}\bm{1}_{m}]_{+}$, where

The expression for $\bm{\alpha}(\bm{\beta})$ is completely symmetrical. While a projection onto the simplex is required for each coordinate, as discussed in §3.3, this can be done in expected linear time. In addition, each coordinate-wise solution can be computed in parallel.

Alternating minimization. Once we know how to compute $\bm{\beta}(\bm{\alpha})$ and $\bm{\alpha}(\bm{\beta})$, there are a number of ways we can build a proper algorithm to solve the smoothed dual. Perhaps the simplest is to alternate between $\bm{\beta}\leftarrow\bm{\beta}(\bm{\alpha})$ and $\bm{\alpha}\leftarrow\bm{\alpha}(\bm{\beta})$. For entropic regularization, this two-block coordinate descent (CD) scheme is known as the Sinkhorn algorithm and was recently popularized in the context of optimal transport by Cuturi (2013). A disadvantage of this approach, however, is that computational effort is spent updating coordinates that may already be near-optimal. To address this issue, we can instead adopt a greedy CD scheme as recently proposed for entropic regularization by Altschuler et al. (2017).

Appendix C Additional experiments

We ran the same experiments as Figure 2 and Figure 3 on one more image pair: “Grafiti” by Jon Ander and “Rainbow Bridge National Monument Utah”, by Bernard Spragg. Both images are in the public domain. The results, presented in Figure 5 and Figure 6 below, confirm the empirical findings described in §6.1 and §6.2. The images are available at https://github.com/mblondel/smooth-ot/tree/master/data.