Arnaud Dessein Qucit

Nicolas Papadakis
email: nicolas.papadakis@math.u-bordeaux.fr

Jean-Luc Rouas
email: jean-luc.rouas@labri.fr

Regularized Optimal Transport and the ROT Mover's Distance

This paper presents a unified framework for smooth convex regularization of discrete optimal transport problems. In this context, the regularized optimal transport turns out to be equivalent to a matrix nearness problem with respect to Bregman divergences. Our framework thus naturally generalizes a previously proposed regularization based on the Boltzmann-Shannon entropy related to the Kullback-Leibler divergence, and solved with the Sinkhorn-Knopp algorithm. We call the regularized optimal transport distance the rot mover's distance in reference to the classical earth mover's distance. By exploiting alternate Bregman projections, we develop the alternate scaling algorithm and non-negative alternate scaling algorithm, to compute efficiently the regularized optimal plans depending on whether the domain of the regularizer lies within the non-negative orthant or not. We further enhance the separable case with a sparse extension to deal with high data dimensions. We also instantiate our framework and discuss the inherent specificities for well-known regularizers and statistical divergences in the machine learning and information geometry communities. Finally, we demonstrate the merits of our methods with experiments using synthetic data to illustrate the effect of different regularizers, penalties and dimensions, as well as real-world data for a pattern recognition application to audio scene classification.

++

to produce a contingency table π that has fixed marginals p, q by Kullback-Leibler projection of ξ onto Π(p, q): π = argmin

They showed that alternate Bregman projections specialize to the Sinkhorn-Knopp algorithm in this context. However, no relationship to OT problems was highlighted.

Contributions and Organization

Our main contribution is to formulate a unified framework for discrete regularized optimal transport (ROT) by considering a large class of smooth convex regularizers. We call the underlying distance the rot mover's distance (RMD) and show that a given ROT problem actually amounts to the minimization of an associated Bregman divergence. This allows the derivation of two schemes that we call the alternate scaling algorithm (ASA) and the non-negative alternate scaling algorithm (NASA), to compute efficiently the regularized optimal plans depending on whether the domain of the regularizer lies within the non-negative orthant or not. These schemes are based on the general form of alternate projections for Bregman divergences. They also exploit the Newton-Raphson method to approximate the projections for separable divergences. The separable case is further enhanced with a sparse extension to deal with high data dimensions. We also instantiate our two generic schemes with widely-used regularizers and statistical divergences.

, where Bregman divergences are known to possess a dually flat structure with a generalized Pythagorean theorem in relation to information projections.

The remainder of this paper is organized as follows. In Section 2, we introduce some necessary preliminaries. In Section 3, we present our theoretical results for a unified framework of ROT problems. We then derive the algorithmic methods for solving ROT problems in Section 4. We also discuss the inherent specificities of ROT problems for classical regularizers and associated divergences in Section 5. In Section 6, we provide experiments to illustrate our methods on synthetic data and real-world audio data in a classification problem. Finally, in Section 7, we draw some conclusions and perspectives for future work.

Introduction

A recurrent problem in statistical machine learning is the choice of a relevant distance measure to compare probability distributions. Various information divergences are famous, among which Euclidean, Mahalanobis, Kullback-Leibler, Itakura-Saito, Hellinger, χ 2 , p (quasi-)norm, total variation, logistic loss function, or more general Csiszár and Bregman divergences and parametric families of such divergences such as α-and β-divergences.

An alternative family of distances between probability distributions can be introduced in the framework of optimal transport (OT). Rather than performing a pointwise comparison of the distributions, the idea is to quantify the minimal effort for moving the probability mass of one distribution to the other, where the transport plan to move the mass is optimized according to a given ground cost. This makes OT distances suitable and robust in certain applications, notably in the field of computer vision where the discrete OT distance, also known as earth mover's distance (EMD), has been popularized to compare histograms of features for pattern recognition tasks [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF].

Despite its appealing theoretical properties, intuitive formulation, and excellent performance in various problems of information retrieval, the computation of the EMD involves solving a linear program whose cost quickly becomes prohibitive with the data dimension. In practice, the best algorithms currently proposed, such as the network simplex [START_REF] Ahuja | Network Flows: Theory, Algorithms and Applications[END_REF], scale at least with a super-cubic complexity. Embeddings of the distributions can be used to approximate the EMD with linear complexity [START_REF] Grauman | Fast contour matching using approximate earth mover's distance[END_REF][START_REF] Indyk | Fast image retrieval via embeddings[END_REF][START_REF] Shirdhonkar | Approximate earth mover's distance in linear time[END_REF], and the network simplex can be modified to run in quadratic time [START_REF] Gudmundsson | Small Manhattan networks and algorithmic applications for the earth mover's distance[END_REF][START_REF] Ling | An efficient earth mover's distance algorithm for robust histogram comparison[END_REF][START_REF] Pele | Fast and robust earth mover's distances[END_REF]. Nevertheless, the distortions inherent to such embeddings [START_REF] Naor | Planar earthmover is not in l 1[END_REF], and the exponential increase of costs incurred by such modifications, make these approaches inapplicable for dimensions higher than four. Instead, multi-scale strategies [START_REF] Oberman | An efficient linear programming method for optimal transportation[END_REF] and shortcut paths [START_REF] Schmitzer | A sparse multi-scale algorithm for dense optimal transport[END_REF] can speed up the estimation of the exact optimal plan. These approaches are yet limited to particular convex costs such as 2 , while other costs such as 1 and truncated or compressed versions are often preferred in practice for an increased robustness to data outliers [START_REF] Pele | A linear time histogram metric for improved SIFT matching[END_REF][START_REF] Pele | Fast and robust earth mover's distances[END_REF][START_REF] Rabin | A statistical approach to the matching of local features[END_REF]. For general applications, a gain in performance can also be obtained with a cost directly learned from labeled data [START_REF] Cuturi | Ground metric learning[END_REF]. The aforementioned accelerated methods that are dedicated to 2 or convex costs are thus not adapted in this context.

On another line of research, the regularization of the transport plan, for example via graph modeling [START_REF] Ferradans | Regularized discrete optimal transport[END_REF], has been considered to deal with noisy data, though this latter approach does not address the computational issue of efficiency for high dimensions. In this continuity, an entropic regularization was shown to admit an efficient algorithm with quadratic complexity that speeds up the computation of solutions by several orders of magnitude, and to improve performance on applications such as handwritten digit recognition [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF]. In addition, a tailored computation can be obtained via convolution for specific ground costs [START_REF] Solomon | Convolutional Wasserstein distances: Efficient optimal transportation on geometric domains[END_REF]. Since the introduction of the entropic regularization, OT has benefited from extensive developments in the machine learning community, with applications such as label propagation [START_REF] Solomon | Wasserstein propagation for semi-supervised learning[END_REF], domain adaptation [START_REF] Courty | Optimal transport for domain adaptation[END_REF], matrix factorization [START_REF] Zen | Simultaneous ground metric learning and matrix factorization with earth mover's distance[END_REF], dictionary learning [START_REF] Rolet | Fast dictionary learning with a smoothed Wasserstein loss[END_REF][START_REF] Schmitz | Wasserstein dictionary learning: Optimal transport-based unsupervised non-linear dictionary learning[END_REF], barycenter computation [START_REF] Cuturi | A smoothed dual approach for variational Wasserstein problems[END_REF], geodesic principal component analysis [START_REF] Bigot | Geodesic PCA in the Wasserstein space[END_REF][START_REF] Cazelles | Log-PCA versus geodesic PCA of histograms in the Wasserstein space[END_REF][START_REF] Seguy | Principal geodesic analysis for probability measures under the optimal transport metric[END_REF], data fitting [START_REF] Frogner | Learning with a Wasserstein loss[END_REF], statistical inference [START_REF] Bernton | Inference in generative models using the Wasserstein distance[END_REF], training of Boltzmann machines [START_REF] Montavon | Wasserstein training of restricted Boltzmann machines[END_REF] and generative adversarial networks [4,[START_REF] Bousquet | From optimal transport to generative modeling: the VEGAN cookbook[END_REF][START_REF] Genevay | GAN and VAE from an optimal transport point of view[END_REF].

With the entropic regularization, the gain in computational time is only important for high dimensions or large levels of regularization. For low regularization, advanced optimization strategies can still be used to obtain a significant speed-up [START_REF] Schmitz | Wasserstein dictionary learning: Optimal transport-based unsupervised non-linear dictionary learning[END_REF][START_REF] Thibault | Overrelaxed sinkhorn-knopp algorithm for regularized optimal transport[END_REF]. It is also a well-known effect that the entropic regularization overspreads the transported mass, which may be undesirable for certain applications as in the case of interpolation purposes. An interesting perspective of these works, however, is that many more regularizers are worth investigating to solve OT problems both efficiently and robustly [START_REF] Blondel | Smooth and sparse optimal transport[END_REF][START_REF] Galichon | Cupid's invisible hand: Social surplus and identification in matching models[END_REF][START_REF] Muzellec | Tsallis regularized optimal transport and ecological inference[END_REF]. This is the idea we address in the present work, focusing on smooth convex regularization.

Notations

For the sake of simplicity, we consider distributions with same dimension d, and thus work with the Euclidean space R d×d of square matrices. It is straightforward, however, to extend all results for a different number of bins m, n by using rectangular matrices in R m×n instead. We denote the null matrix of R d×d by 0, and the matrix full of ones by 1. The Frobenius inner product between two matrices π, ξ ∈ R d×d is defined by:

π, ξ = d i=1 d j=1 π ij ξ ij .
(1)

When the intended meaning is clear from the context, we also write 0 for the null vector of R d , and 1 for the vector full of ones. The notation • represents the transposition operator for matrices or vectors. The probability simplex of R d is defined as follows:

Σ d = {p ∈ R d + : p 1 = 1} . ( 2
)
The operator diag(v) transforms a vector v ∈ R d into a diagonal matrix π ∈ R d×d such that π ii = v i , for all 1 ≤ i ≤ d. The operator vec(π) transforms a matrix π ∈ R d×d into a vector x ∈ R d 2 such that x i+(j-1)d = π ij , for all 1 ≤ i, j ≤ d. The operator sgn(x) for x ∈ R returns -1, 0, +1, if x is negative, null, positive, respectively. Functions of a real variable, such as the absolute value, sign, exponential or power functions, are considered element-wise when applied to matrices. The max operator and inequalities between matrices should also be interpreted element-wise. Matrix divisions are similarly considered element-wise, whereas element-wise matrix multiplications, also known as Hadamard or Schur products, are denoted by to remove any ambiguity with standard matrix multiplications. Lastly, addition or subtraction of a scalar and a matrix should be understood element-wise by replicating the scalar.

Background and Related Work

Given two probability vectors p, q ∈ Σ d , and a cost matrix γ ∈ R d×d + whose coefficients γ ij represent the cost of moving the mass from bin p i to q j , the total cost of a given transport plan, or coupling, π ∈ Π(p, q) can be quantified as π, γ . An optimal cost is then obtained by solving a linear program: d γ (p, q) = min

π∈Π(p,q)

π, γ ,

with the transport polytope of p and q, also known as the polytope of couplings between p and q, defined as the following polyhedron:

Π(p, q) = {π ∈ R d×d + : π1 = p, π 1 = q} . ( 4
)
The EMD associated to the cost matrix γ is given by d γ and is a true distance metric on the probability simplex Σ d whenever γ is itself a distance matrix. In general, the optimal plans, or earth mover's plans, have at most 2d -1 nonzero entries, and consist either of a single vertex or of a whole facet of the transport polytope. One of the earth mover's plans can be obtained with the network simplex [START_REF] Ahuja | Network Flows: Theory, Algorithms and Applications[END_REF] among other approaches. For a general cost matrix γ, the complexity of solving an OT problem scales at least in O(d 3 log d) for the best algorithms currently proposed, including the network simplex, and turns out to be super-cubic in practice as well. [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] proposed a new family of OT distances, called Sinkhorn distances, from the perspective of maximum entropy. The idea is to smooth the original problem with a strictly convex regularization via the Boltzmann-Shannon entropy. The primal problem involves the entropic regularization as an additional constraint: d γ,α (p, q) = min π∈Πα(p,q) π, γ , [START_REF] Bauschke | Dykstra's algorithm with Bregman projections: A convergence proof[END_REF] with the regularized transport polytope defined as follows:

Π α (p, q) = {π ∈ Π(p, q) : E(π) ≤ E(pq ) + α} , [START_REF] Benamou | Iterative Bregman projections for regularized transportation problems[END_REF] where α ≥ 0 is a regularization term and E is minus the Boltzmann-Shannon entropy as defined in [START_REF] Indyk | Fast image retrieval via embeddings[END_REF]. It is also straightforward to prove that we have: Π α (p, q) = {π ∈ Π(p, q) : K(π 1) ≤ K(pq 1) + α} , [START_REF] Bernton | Inference in generative models using the Wasserstein distance[END_REF] Π α (p, q) = {π ∈ Π(p, q) : K(π pq ) ≤ α} . [START_REF] Bigot | Geodesic PCA in the Wasserstein space[END_REF] where K is the Kullback-Leibler divergence as defined in [START_REF] Idel | A review of matrix scaling and Sinkhorn's normal form for matrices and positive maps[END_REF]. This enforces the solution to have sufficient entropy, or equivalently small enough mutual information, by constraining it to the Kullback-Leibler ball of radius K(pq 1) + α, respectively α, and center the matrix 1 ∈ R d×d ++ , respectively the transport plan pq ∈ R d×d ++ , which have maximum entropy. The dual problem exploits a Lagrange multiplier to relax the entropic regularization as a penalty:

d γ,λ (p, q) = π λ , γ , (9)
with the regularized optimal plan π λ defined as follows:

π λ = argmin π∈Π(p,q)
π, γ + λE(π) , [START_REF] Bonneel | Displacement interpolation using lagrangian mass transport[END_REF] where λ > 0 is a regularization term. The problem can then be solved empirically in quadratic complexity with linear convergence using the Sinkhorn-Knopp algorithm [START_REF] Sinkhorn | Concerning nonnegative matrices and doubly stochastic matrices[END_REF] based on iterative matrix scaling, where rows and columns are rescaled in turn so that they respectively sum up to p and q until convergence. Finally, it is easy to prove that we have:

π λ = argmin π∈Π(p,q)
π, γ + λK(π 1) , [START_REF] Bousquet | From optimal transport to generative modeling: the VEGAN cookbook[END_REF] π λ = argmin π∈Π(p,q) π, γ + λK(π pq ) . [START_REF] Cazelles | Log-PCA versus geodesic PCA of histograms in the Wasserstein space[END_REF] This again shows that the regularization enforces the solution to have sufficient entropy, or equivalently small enough mutual information, by shrinking it toward the matrix 1 and the joint distribution pq which have maximum entropy. [START_REF] Benamou | Iterative Bregman projections for regularized transportation problems[END_REF] revisited the entropic regularization in a geometrical framework with iterative information projections. They showed that computing a Sinkhorn distance in dual form actually amounts to the minimization of a Kullback-Leibler divergence:

π λ = argmin π∈Π(p,q)
K(π exp(-γ/λ)) . [START_REF] Cominetti | Asymptotic analysis of the exponential penalty trajectory in linear programming[END_REF] Precisely, this amounts to computing the Kullback-Leibler projection of exp(-γ/λ) ∈ R d×d ++ onto the transport polytope Π(p, q). In this context, the Sinkhorn-Knopp algorithm turns out to be a special instance of Bregman projection onto the intersection of convex sets via alternate projections. Specifically, we see Π(p, q) as the intersection of the non-negative orthant with two affine subspaces containing all matrices with rows and columns summing to p and q respectively, and we alternate projection on these two subspaces according to the Kullback-Leibler divergence until convergence.

[29] further studied this equivalence in the wider context of iterative proportional fitting. He notably showed that the Sinkhorn-Knopp and alternate Bregman projections can be extended to account for infinite entries in the cost matrix γ, and thus null entries in the regularized optimal plan. Hence, it is possible to develop a sparse version of the entropic regularization to OT problems. This becomes interesting to store the d × d matrix variables and perform the required computations when the data dimension gets large. [START_REF] Dhillon | Matrix nearness problems with Bregman divergences[END_REF] had already enlightened such an equivalence in the field of matrix analysis. They actually considered the estimation of contingency tables with fixed marginals as a matrix nearness problem based on the Kullback-Leibler divergence. In more detail, they use a rough estimate ξ ∈ R d×d

Theoretical Preliminaries

In this section, we introduce the required preliminaries to our framework. We begin with elements of convex analysis (Section 2.1) and of Bregman geometry (Section 2.2). We proceed with theoretical results for convergence of alternate Bregman projections (Section 2.3) and of the Newton-Raphson method (Section 2.4).

Convex Analysis

Let E be a Euclidean space with inner product •, • and induced norm • . The boundary, interior and relative interior of a subset X ⊆ E are respectively denoted by bd(X ), int(X ), and ri(X ), where we recall that for a convex set C, we have:

ri(C) = {x ∈ C : ∀y ∈ C, ∃λ > 1, λx + (1 -λ)y ∈ C} . ( 15
)
In convex analysis, scalar functions are defined over the whole space E and take values in the extended real number line R ∪ {-∞, +∞}. The effective domain, or simply domain, of a function f is then defined as the set:

dom f = {x ∈ E : f (x) < +∞} . ( 16
)
A convex function f is proper if f (x) < +∞ for at least one x ∈ E and f (x) > -∞ for all x ∈ E, and it is closed if its lower level sets {x ∈ E : f (x) ≤ α} are closed for all α ∈ R. If dom f is closed, then f is closed, and a proper convex function is closed if and only if it is lower semi-continuous. Moreover, a closed function f is continuous relative to any simplex, polytope or polyhedral subset in dom f . It is also well-known that a convex function f is always continuous in the relative interior ri(dom f ) of its domain.

A function f is essentially smooth if it is differentiable on int(dom f ) = ∅ and verifies lim k→+∞ ∇f (x k ) = +∞ for any sequence (x k ) k∈N from int(dom f ) that converges to a point x ∈ bd(dom f ). A function f is of Legendre type if it is a closed proper convex function that is also essentially smooth and strictly convex on int(dom f ).

The Fenchel conjugate f of a function f is defined for all y ∈ E as follows:

f (y) = sup x∈int(dom f ) x, y -f (x) . ( 17
)
The Fenchel conjugate f is always a closed convex function. Moreover, if f is a closed convex function, then (f ) = f , and f is of Legendre type if and only if f is of Legendre type. In this latter case, the gradient mapping ∇f is a homeomorphism between int(dom f ) and int(dom f ), with inverse mapping (∇f ) -1 = ∇f , which guarantees the existence of dual coordinate systems x(y) = ∇f (y) and y(x) = ∇f (x) on int(dom f ) and int(dom f ).

Finally, we say that a function f is cofinite if it verifies:

lim λ→+∞ f (λx)/λ = +∞ , ( 18
)
for all nonzero x ∈ E. Intuitively, it means that f grows super-linearly in every direction. In particular, a closed proper convex function is cofinite if and only if dom f = E.

Bregman Geometry

Let φ be a convex function on E that is differentiable on int(dom φ) = ∅. The Bregman divergence generated by φ is defined as follows:

B φ (x y) = φ(x) -φ(y) -x -y, ∇φ(y) , ( 19
)
for all x ∈ dom φ and y ∈ int(dom φ). We have B φ (x y) ≥ 0 for any x ∈ dom φ and y ∈ int(dom φ). If in addition φ is strictly convex on int(dom φ), then B φ (x y) = 0 if and only if x = y. Bregman divergences are also always convex in the first argument, and are invariant under adding an arbitrary affine term to their generator.

Bregman divergences are not symmetric and do not verify the triangle inequality in general, and thus are not necessarily distances in the strict sense. However, they still enjoy some nice geometrical properties that somehow generalize the Euclidean geometry. In particular, they verify a four-point identity similar to a parallelogram law:

B φ (x y) + B φ (x y ) = B φ (x y) + B φ (x y ) -x -x , ∇φ(y) -∇φ(y ) , ( 20
)
for all x, x ∈ dom φ and y, y ∈ int(dom φ). A special instance of this relation gives rise to a three-point property similar to a triangle law of cosines:

B φ (x y) = B φ (x y ) + B φ (y y) -x -y , ∇φ(y) -∇φ(y ) , ( 21
)
for all x ∈ dom φ and y, y ∈ int(dom φ). Suppose now that φ is of Legendre type, and let C ⊆ E be a closed convex set such that C ∩ int(dom φ) = ∅. Then, for any point y ∈ int(dom φ), the following problem:

P C (y) = argmin x∈C B φ (x y) , ( 22
)
has a unique solution, then called the Bregman projection of y onto C. This solution actually belongs to C ∩ int(dom φ), and is also characterized as the unique point y ∈ C ∩ int(dom φ) that verifies the variational relation:

x -y , ∇φ(y) -∇φ(y ) ≤ 0 , (23)
for all x ∈ C∩dom φ. This characterization is equivalent to a well-known generalized Pythagorean theorem for Bregman divergences, which states that the Bregman projection of y onto C is the unique point y ∈ C ∩ int(dom φ) that verifies the following inequality:

B φ (x y) ≥ B φ (x y ) + B φ (y y) , ( 24
)
for all x ∈ C ∩ dom φ. When C is further an affine subspace, or more generally when the Bregman projection further belongs to ri(C), the scalar product actually vanishes:

xy , ∇φ(y) -∇φ(y ) = 0 , [START_REF] Gudmundsson | Small Manhattan networks and algorithmic applications for the earth mover's distance[END_REF] leading to an equality in the generalized Pythagorean theorem:

B φ (x y) = B φ (x y ) + B φ (y y) . ( 26
)
A famous example of Bregman divergence is the Kullback-Leibler divergence, defined for matrices π ∈ R d×d + and ξ ∈ R d×d ++ as follows:

K(π ξ) = d i=1 d j=1 π ij log π ij ξ ij -π ij + ξ ij . ( 27
)
This divergence is generated by a function of Legendre type for π ∈ R d×d

+

given by minus the Boltzmann-Shannon entropy:

E(π) = K(π 1) = d i=1 d j=1 (π ij log(π ij ) -π ij + 1) , ( 28
)
with the convention 0 log(0) = 0. Another well-known example is the Itakura-Saito divergence, defined for matrices π, ξ ∈ R d×d ++ as follows:

I(π ξ) = d i=1 d j=1 π ij ξ ij -log π ij ξ ij -1 . ( 29
)
This divergence is generated by a function of Legendre type for π ∈ R d×d ++ given by minus the Burg entropy:

F (π) = d i=1 d j=1 (π ij -log π ij -1) . ( 30
)
On the one hand, these examples belong to a particular type of so-called separable Bregman divergences between matrices on R d×d , that can be seen as the aggregation of element-wise Bregman divergences between scalars on R:

B φ (π ξ) = d i=1 d j=1 B φij (π ij ξ ij ) , ( 31
)
φ(π) = d i=1 d j=1 φ ij (π ij ) . ( 32
)
Often, all element-wise generators φ ij are chosen equal, and are thus simply written as φ with a slight abuse of notation. Other examples of such divergences are discussed in Section 5, and include the logistic loss function generated by minus the Fermi-Dirac entropy, or the squared Euclidean distance generated by the Euclidean norm.

On the other hand, a classical example of non-separable Bregman divergence is half the squared Mahalanobis distance, defined for matrices π, ξ ∈ R d×d as follows:

M (π ξ) = 1 2 vec(π -ξ) P vec(π -ξ) , ( 33
)
for a positive-definite matrix P ∈ R d 2 ×d 2 . This divergence is generated by a function of Legendre type for π ∈ R d×d given by a quadratic form:

Q(π) = 1 2 vec(π) P vec(π) . ( 34
)
This example is also discussed in Section 5.

Alternate Bregman Projections

Let φ be a function of Legendre type with Fenchel conjugate φ = ψ. In general, computing Bregman projections onto an arbitrary closed convex set C ⊆ E such that C ∩ int(dom φ) = ∅ is nontrivial. Sometimes, it is possible to decompose C into the intersection of finitely many closed convex sets:

C = s l=1 C l , ( 35
)
where the individual Bregman projections onto the respective sets C 1 , . . . , C s are easier to compute. It is then possible to obtain the Bregman projection onto C by alternate projections onto C 1 , . . . , C s according to Dykstra's algorithm.

In more detail, let σ : N → {1, . . . , s} be a control mapping that determines the sequence of subsets onto which we project. For a given point x 0 ∈ C ∩ int(dom φ), the Bregman projection P C (x 0 ) of x 0 onto C can be approximated with Dykstra's algorithm by iterating the following updates:

x k+1 ← P C σ(k) (∇ψ(∇φ(x k ) + y σ(k) )) , ( 36
)
where the correction terms y 1 , . . . , y s for the respective subsets are initialized with the null element of E, and are updated after projection as follows:

y σ(k) ← y σ(k) + ∇φ(x k ) -∇φ(x k+1 ) . ( 37
)
Under some technical conditions, the sequence of updates (x k ) k∈N then converges in norm to P C (x 0 ) with a linear rate. Several sets of such conditions have been studied, notably by [START_REF] Tseng | Dual coordinate ascent methods for non-strictly convex minimization[END_REF], [START_REF] Bauschke | Dykstra's algorithm with Bregman projections: A convergence proof[END_REF], [START_REF] Dhillon | Matrix nearness problems with Bregman divergences[END_REF].

We here use the conditions proposed by [START_REF] Dhillon | Matrix nearness problems with Bregman divergences[END_REF], which reveal to be the less restrictive ones in our framework. Specifically, the convergence of Dykstra's algorithm is guaranteed as soon as the function φ is cofinite, the constraint qualification ri(C 1 ) ∩ • • • ∩ ri(C s ) ∩ int(dom φ) = ∅ holds, and the control mapping σ is essentially cyclic, that is, there exists a number t ∈ N such that σ takes each output value at least once during any t consecutive input values. If a given C l is a polyhedral set, then the relative interior can be dropped from the constraint qualification. Hence, when all subsets C l are polyhedral, the constraint qualification simply reduces to C ∩int(dom φ) = ∅, which is already enforced for the definition of Bregman projections.

Finally, if all subsets C l are further affine, then we can relax other assumptions. Notably, we do not require φ to be cofinite [START_REF] Dessein | Parameter estimation in finite mixture models by regularized optimal transport: A unified framework for hard and soft clustering[END_REF], or equivalently dom ψ = E, but only dom ψ to be open. The control mapping need not be essentially cyclic anymore, as long as it takes each output value an infinite number of times. More importantly, we can completely drop the correction terms from the updates, leading to a simpler technique known as projections onto convex sets (POCS):

x k+1 ← P C σ(k) (x k ) . ( 38
)

Newton-Raphson Method

Let f be a continuously differentiable scalar function on an open interval I ⊆ R. Assume f is increasing on a non-empty closed interval [x -, x + ] ⊂ I, and write y -= f (x -) and y + = f (x + ). Then, for any y ∈ [y -, y + ], the equation f (x) = y has at least one solution x ∈ [x -, x + ]. Such a solution can be approximated by iterative updates according to the Newton-Raphson method:

x ← max x -, min x + , x - f (x) -y f (x) , ( 39
)
where the fraction takes infinite values when f (x) = 0 and f (x) = y, and a null value by convention when f (x) = 0 and f (x) = y. It is well-known that the Newton-Raphson method converges to a solution x as soon as x is initialized sufficiently close to x . Convergence is then quadratic provided that f (x ) = 0. However, this local convergence has little importance in practice because it is hard to quantify the required proximity to the solution.

[52] elucidated results on global convergence of the Newton-Raphson method. He proved a necessary and sufficient condition of convergence for an arbitrary value y ∈ [y -, y + ] and from any starting point

x ∈ [x -, x + ]. This condition is that for any a, b ∈ [x -, x + ], f (b) > f (a) implies: f (a) + f (b) > f (b) -f (a) b -a . ( 40
) (A) Affine constraints (B) Polyhedral constraints (A1) φ is of Legendre type. (B1) φ is of Legendre type. (A2) (0, 1) d×d ⊆ dom φ. (B2) (0, 1) d×d ⊆ dom φ. (A3) dom φ ⊆ R d×d + . (B3) dom φ R d×d + . (A4) dom ψ is open. (B4) dom ψ = R d×d . (A5) R d×d - ⊂ dom ψ.
Table 1: Set of assumptions for the considered regularizers φ.

In particular, a sufficient condition is that the underlying function f is an increasing convex or increasing concave function on [x -, x + ], or can be decomposed as the sum of such functions. In addition, if f satisfies the necessary and sufficient condition and is strictly increasing with f (x) > 0 for all x ∈ [x -, x + ], then initializing with a boundary point x -= x or x + = x ensures that the entire sequence of updates is interior to (x -, x + ), so that we can actually drop the min and max truncation operators in the updates:

x ← x - f (x) -y f (x) . ( 41
)

Mathematical Formulation

In this section, we develop a unified framework to define ROT problems. We start by drawing some technical assumptions for our generalized framework to hold (Section 3.1). We then formulate primal ROT problems and study their properties (Section 3.2). We also formulate dual ROT problems and discuss their properties in relation to primal ones (Section 3.3). Finally, we provide some geometrical insights to summarize our developments in the light of information geometry (Section 3.4).

Technical Assumptions

Some mild technical assumptions are required on the convex regularizer φ and its Fenchel conjugate ψ = φ * for the proposed framework to hold. Some assumptions relate to required conditions for the definition of Bregman projections and convergence of the algorithms, while others are more specific to ROT problems. In our framework, we also need to distinguish between two situations where the underlying closed convex set can be described as the intersection of either affine subspaces or polyhedral subsets. The two sets of assumptions (A) and (B) are summarized in Table 1.

For the first assumptions (A1) and (B1), we recall that a closed proper convex function is of Legendre type if and only if it is essentially smooth and strictly convex on the interior of its domain (Section 2.1). This is required for the definition of Bregman projections (Section 2.2). In addition, it guarantees the existence of dual coordinate systems on int(dom φ) and int(dom ψ) via the homeomorphism ∇φ = ∇ψ -1 :

π(θ) = ∇ψ(θ) , ( 42
)
θ(π) = ∇φ(π) . ( 43
)
With a slight abuse of notation, we omit the reparameterization to simply denote corresponding primal and dual parameters by π and θ.

The second assumptions (A2) and (B2) imply that ri(Π(p, q)) ⊂ dom φ and ensure the constraint qualification Π(p, q) ∩ int(dom φ) = ∅ for Bregman projection onto the transport polytope, independently of the input distributions p, q as long as they do not have null or unit entries. We assume hereafter that this implicitly holds, and discuss in the practical considerations (Section 4.6) how our methods actually generalize to deal explicitly with null or unit entries in the input distributions.

The third assumptions (A3) and (B3) separate between two cases depending on whether dom φ lies within the non-negative orthant or not for the alternate Bregman projections (Section 2.3). In the former case, non-negativity is already ensured by the domain of the regularizer, so that the underlying closed convex set is made of two affine subspaces for the row and column sum constraints, and the POCS method can be considered. The fourth assumption (A4) thus requires that dom ψ be open for convergence of this algorithm. In the latter case, there is one additional polyhedral subset for the non-negative constraints and Dykstra's algorithm should be used. The fourth assumption (B4) hence further requires that dom ψ = R d×d , or equivalently that φ be cofinite [START_REF] Dessein | Parameter estimation in finite mixture models by regularized optimal transport: A unified framework for hard and soft clustering[END_REF], for convergence. In both cases, we remark that we necessarily have dom ψ = dom ∇ψ.

The fifth assumption (A5) in the affine constraints ensures that -γ/λ belongs to dom ∇ψ for definition of ROT problems, independently of the non-negative cost matrix γ and positive regularization term λ. Notice that this is already guaranteed by the fourth assumption in the polyhedral constraints. We also show in the sparse extension (Section 4.5) how to deal with infinite entries in the cost matrix γ for separable regularizers, so as to enforce null entries in the regularized optimal plan.

On the one hand, some common regularizers under assumptions (A) are the Boltzmann-Shannon entropy associated to the Kullback-Leibler divergence, the Burg entropy associated to the Itakura-Saito divergence, and the Fermi-Dirac entropy associated to the logistic loss function. To solve the underlying ROT problems, we employ our method called ASA based on the POCS technique, where alternate Bregman projections onto the two affine subspaces for the row and column sum constraints are considered (Section 4.3). On the other hand, examples under assumptions (B) include the Euclidean norm associated to the Euclidean distance, and the quadratic form associated to the Mahalanobis distance. For these ROT problems, we use our second method called NASA based on Dykstra's algorithm, where correction terms and a further Bregman projection onto the polyhedral non-negative orthant are needed (Section 4.4).

Primal Problem

We start our primal formulation with the following lemmas and definition for the RMD. Lemma 1. The regularizer φ attains its global minimum uniquely at ξ = ∇ψ(0).

Proof. Using the assumptions (A4) and (A5), respectively (B4), we have that 0 ∈ dom ψ = int(dom ψ). Thus, there exists a unique ξ ∈ int(dom φ) such that ∇φ(ξ ) = 0, or equivalently ξ = ∇ψ(0), via the homeomorphism ∇φ = ∇ψ -1 ensured by assumption (A1), respectively (B1). Hence, φ attains its global minimum uniquely at ξ by strict convexity on int(dom φ). Lemma 2. The restriction of the regularizer φ to the transport polytope Π(p, q) attains its global minimum uniquely at the Bregman projection π of ξ onto Π(p, q). Proof. Using the assumption (A2), respectively (B2), we have that Π(p, q) ∩ int(dom φ) = ∅. Since ξ ∈ int(dom φ) and Π(p, q) is a closed convex set, the Bregman projection π of ξ onto Π(p, q) according to the function φ of Legendre type is well-defined. Moreover, it is characterized by the variational relation [START_REF] Genevay | GAN and VAE from an optimal transport point of view[END_REF] as follows:

π -π , ∇φ(π ) ≥ 0 , ( 44
)
for all π ∈ Π(p, q) ∩ dom φ. We also have B φ (π π ) > 0 when π = π by strict convexity of φ on int(dom φ). As a result, we have:

φ(π) -φ(π ) > π -π , ∇φ(π ) . ( 45
)
Combining the two inequalities, we obtain φ(π) > φ(π ) and the restriction of φ to Π(p, q) attains its global minimum uniquely at π . Lemma 3. The restriction of the cost •, γ to the regularized transport polytope:

Π α,φ (p, q) = {π ∈ Π(p, q) : φ(π) ≤ φ(π ) + α} , ( 46
)
where α ≥ 0, attains its global minimum.

Proof. The regularized transport polytope is the intersection of the compact set Π(p, q) with a lower level set of φ which is also closed since φ is closed. Hence, Π α,φ (p, q) is compact and the restriction of •, γ to Π α,φ (p, q) attains its global minimum by continuity on a compact set.

Definition 1. The primal rot mover's distance is the quantity defined as:

d γ,α,φ (p, q) = min π∈Π α,φ (p,q) π, γ . ( 47
)
A minimizer π α is then called a primal rot mover's plan.

Remark 1. For the sake of notation, we omit the dependence on p, q, γ, φ in the index of primal rot mover's plans π α .

The regularization enforces the associated minimizers to have small enough Bregman information φ(π α ) ≤ φ(π ) + α compared to the minimal one φ(π ) for transport plans. We also have a geometrical interpretation where the solutions are constrained to a Bregman ball whose center ξ is the matrix with minimal Bregman information. Proposition 4. The regularized transport polytope is the intersection of the transport polytope with the Bregman ball of radius B φ (π ξ ) + α and center ξ :

Π α,φ (p, q) = {π ∈ Π(p, q) : B φ (π ξ ) ≤ B φ (π ξ ) + α} . ( 48
)
Proof. Expanding the Bregman divergences from their definition [START_REF] Dhillon | Matrix nearness problems with Bregman divergences[END_REF], we obtain:

B φ (π ξ ) = φ(π) -φ(ξ ) -π -ξ , ∇φ(ξ ) , ( 49
)
B φ (π ξ ) = φ(π ) -φ(ξ ) -π -ξ , ∇φ(ξ ) . ( 50
)
Since ∇φ(ξ ) = 0, the last terms with scalar products vanish, leading to:

φ(π) -φ(π ) = B φ (π ξ ) -B φ (π ξ ) . ( 51
)
Therefore, in the definition [START_REF] Sha | Multiplicative updates for nonnegative quadratic programming[END_REF] of Π α,φ (p, q), we have φ(π) ≤ φ(π ) + α if and only if π is in the Bregman ball of radius B φ (π ξ ) + α and center ξ .

Under some additional conditions, this geometrical interpretation still holds with a Bregman ball whose center π has minimal Bregman information for transport plans. Proposition 5. If π ∈ ri(Π(p, q)), then the regularized transport polytope is the intersection of the transport polytope with the Bregman ball of radius α and center π :

Π α,φ (p, q) = {π ∈ Π(p, q) : B φ (π π ) ≤ α} . ( 52
)
Proof. Since π ∈ ri(Π(p, q)), there is equality in the generalized Pythagorean theorem [START_REF] Haasdonk | Feature space interpretation of SVMs with indefinite kernels[END_REF]:

B φ (π ξ ) = B φ (π π ) + B φ (π ξ ) . ( 53
)
The regularized transport polytope as seen from ( 48) is then the intersection of the transport polytope Π(p, q) with the Bregman ball of radius α and center π .

Remark 2. The proposition also holds trivially when the global minimum is attained on the transport polytope, that is, when ξ = π .

Corollary 6. Under assumptions (A), the regularized transport polytope is the intersection of the transport polytope with the Bregman ball of radius α and center π :

Π α,φ (p, q) = {π ∈ Π(p, q) : B φ (π π ) ≤ α} . ( 54
)
Proof. This is a result of π ∈ Π(p, q)∩int(dom φ) = ri(Π(p, q)) when dom φ ⊆ R d×d + . Indeed, we then have ri(Π(p, q)) ⊂ Π(p, q) and ri(Π(p, q)) ⊂ int(dom φ), so that ri(Π(p, q)) ⊆ Π(p, q) ∩ int(dom φ). Conversely, let π ∈ Π(p, q) ∩ int(dom φ) so that π ∈ R d×d ++ . Then, for a given π ∈ Π(p, q), let us pose π λ = λπ + (1 -λ)π for λ > 1. We easily have π λ 1 = p and π λ 1 = q. Moreover, since all entries of π are positive and that of π are non-negative, we can always choose a given λ sufficiently close to 1 such that π λ ∈ R d×d + . We then have π λ ∈ Π(p, q) so that π ∈ ri(Π(p, q)) as characterized by [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF], and thus Π(p, q) ∩ int(dom φ) ⊆ ri(Π(p, q)).

Remark 3. Under assumptions (B), the Bregman projection π does not necessarily lie within ri(Π(p, q)). Hence, the geometrical interpretation in terms of a Bregman ball might break down, although the solutions are still constrained to have a small enough Bregman information above that of π .

Although Sinkhorn distances verify the triangular inequality when γ is a distance matrix, thanks to specific chain rules and information inequalities for the Bolzmann-Shannon entropy and Kullback-Leibler divergence, it is not necessarily the case for the RMD with other regularizations, even for separable regularizers. Hence, the RMD does not provide a true distance metric on Σ d in general even if γ is a distance matrix. Nonetheless, the RMD is symmetric as soon as φ is invariant by transposition, which holds for separable regularizers φ ij = φ, and γ is symmetric. We now study some properties of the RMD that hold for general regularizers. d γ,α,φ (p,q) is a decreasing convex and continuous function of α.

Property 1. The primal rot mover's distance

Proof. The fact that it is decreasing is a direct consequence of the regularized transport polytope Π α,φ (p, q) growing with α. The convexity can be proved as follows. Let α 0 , α 1 ≥ 0, and 0 < λ < 1. We pose α λ = (1 -λ)α 0 + λα 1 ≥ 0. We also choose arbitrary rot mover's plans π α0 , π α1 , π α λ . We finally pose π λ = (1 -λ)π α0 + λπ α1 . By convexity of φ, we have:

φ(π λ ) ≤ (1 -λ)φ(π α0 ) + λφ(π α1 ) (55) ≤ (1 -λ)(α 0 + φ(π )) + λ(α 1 + φ(π )) (56) = α λ + φ(π ) . ( 57
)
Hence, π λ ∈ Π α λ ,φ (p, q), and by construction we have π α λ , γ ≤ π λ , γ , or equivalently:

π α λ , γ ≤ (1 -λ) π α0 , γ + λ π α1 , γ . ( 58
)
The continuity for α > 0 is a direct consequence of convexity for α > 0, since a convex function is always continuous on the relative interior of its domain. Lastly, the continuity at α = 0 can be seen as follows. Let (α k ) k∈N be a sequence of positive numbers that converges to 0. We choose arbitrary rot mover's plans (π α k ) k∈N . By compactness of Π(p, q), we can extract a subsequence of rot mover's plans that converges in norm to a point π ∈ Π(p, q). For the sake of simplicity, we do not relabel this subsequence. By construction, we have φ(π

) ≤ φ(π α k ) ≤ φ(π ) + α k ,
and φ(π α k ) converges to φ(π ). By lower semi-continuity of φ, we thus have φ(π ) ≤ φ(π ).

Since the global minimum of φ on Π(p, q) is attained uniquely at π , we must have π = π , and the original sequence also converges in norm to π . By continuity of the total cost •, γ on R d×d , π α k , γ converges to π , γ . Hence, the limit of the RMD when α tends to 0 from above is π , γ , which equals the RMD for α = 0 as shown in the next property.

Property 2. When α = 0, the primal rot mover's distance reduces to:

d γ,0,φ (p, q) = π , γ , ( 59
)
and the unique primal rot mover's plan is the transport plan with minimal Bregman information:

π 0 = π . ( 60
)
Proof. Since π is the unique global minimizer of φ on Π(p, q), the regularized transport polytope reduces to the singleton Π 0,φ (p, q) = {π ∈ Π(p, q) : φ(π) ≤ φ(π )} = {π }. The property follows immediately.

Property 3. When α tends to +∞, the primal rot mover's distance converges to the earth mover's distance:

lim α→+∞ d γ,α,φ (p, q) = d γ (p, q) . ( 61
)
Proof. Let π ∈ Π(p, q) be an earth mover's plan so that d γ (p, q) = π , γ . By continuity of the total cost •, γ on R d×d , we have that for all > 0, there exists an open neighborhood of π such that π, γ ≤ π , γ + for any transport plan π within this neighborhood. We can always choose a transport plan such that π ∈ ri(Π(p, q)). Since ri(Π(p, q)) ⊂ dom φ, φ(π) is finite and we can fix α = φ(π) -φ(π ) ≥ 0. Hence, π ∈ Π α,φ (p, q) for any α ≥ α , and we have

d γ (p, q) ≤ d γ,α,φ (p, q) ≤ π, γ ≤ d γ (p, q) + . Property 4. If [0, 1)
d×d ⊆ dom φ, then there exists a minimal α ≥ 0 such that for all α ≥ α , the primal rot mover's distance reduces to the earth mover's distance:

d γ,α,φ (p, q) = d γ (p, q) . ( 62
)
Proof. The extra condition guarantees that Π(p, q) ⊂ dom φ, and thus that φ is bounded on the closed set Π(p, q). The property is then a direct consequence of Π α,φ (p, q) = Π(p, q) for α large enough.

Property 5. If [0, 1) d×d ⊆ dom φ and φ is strictly convex on [0, 1) d×d , then the unique primal rot mover's plan for α = α is the earth mover's plan π 0 with minimal Bregman information:

π α = π 0 . ( 63
)
Proof. First, we recall that the set of earth mover's plans π is either a single vertex or a whole facet of Π(p, q). Hence, it forms a closed convex subset in Π(p, q), and there is a unique earth mover's plan π 0 with minimal Bregman information by strict convexity of φ on this subset. Second, it is trivial that all primal rot mover's plan π α must be earth mover's plans. If there is a single vertex as earth mover's plan, then the property follows immediately. Otherwise, we can see the property geometrically as follows. The whole facet of earth mover's plans is orthogonal to γ. Nevertheless, by strict convexity of φ on [0, 1) d×d , the facet must be tangent to Π α ,φ (p, q) at the unique earth mover's plan π 0 with minimal Bregman information φ(π 0 ) = φ(π ) + α , and π 0 is also the rot mover's plan π α . Another way to prove the property more formally is as follows. Suppose that a primal rot mover's plan π is not the earth mover's plan with minimal Bregman information. We thus have φ(π 0 ) < φ(π ) ≤ φ(π ) + α . We can then choose a smaller α such that φ(π 0 ) ≤ φ(π ) + α and the RMD still equals the EMD for this smaller value, and actually all values in between by monotonicity. This leads to a contradiction and π 0 must be the earth mover's plan with minimal Bregman information.

Remark 4. When α > α , the regularized transport polytope might grow to include several earth mover's plans with different Bregman information, which are then all minimizers for the RMD. When we do not have strict convexity outside (0, 1)

d×d , there might also be multiple earth mover's plans with minimal Bregman information.

If [0, 1)

d×d ⊆ dom φ, then it is easy to check that the strict convexity of φ on [0, 1) d×d is always verified when φ is separable under assumptions (A) or (B), or when [0, 1) d×d ⊂ int(dom φ) under assumptions (B). This holds for almost all typical regularizers, notably for all regularizers considered in this paper except from minus the Burg entropy as defined in [START_REF] Lin | A study on sigmoid kernels for SVM and the training of non-PSD kernels by SMO-type methods[END_REF] and associated to the Itakura-Saito divergence in [START_REF] Kurras | Symmetric iterative proportional fitting[END_REF]. For this latter regularizer, the solutions for an increasing α all lie within ri(Π(p, q)), and the RMD never reaches the EMD. In such cases where the minimal α does not exist, we can use the convention α = +∞ since the RMD always converges to the EMD in the limit when α tends to +∞. We can then prove that there is a unique rot mover's plan π α as long as 0 < α < α , which can be seen informally as follows. The solutions geometrically lie at the intersection of Π α,φ (p, q) and of a supporting hyperplane with normal γ. By strict convexity of φ on ri(Π(p, q)), this intersection is a singleton inside the polytope. When the intersection reaches a facet, the only facet that can coincide locally with the hyperplane is the one that contains the earth mover's plans. Hence, we also have a singleton on the boundary of the polytope before reaching an earth mover's plan. We formally prove this uniqueness result next by exploiting duality.

Dual Problem

We now present the following two lemmas before defining our dual formulation for the RMD.

Lemma 7. The regularized cost •, γ + λφ(•)

, where λ > 0, attains its global minimum uniquely at ξ = ∇ψ(-γ/λ).

Proof. The regularized cost is convex with same domain as φ, and is strictly convex on int(dom φ). Thus, it attains its global minimum at a unique point ξ ∈ int(dom φ) if and only if γ + λ∇φ(ξ) = 0, or equivalently ∇φ(ξ) = -γ/λ. By assumptions (A4) and (A5), respectively (B4), -γ/λ ∈ dom ∇ψ, so that the global minimum is attained uniquely at ξ = ∇ψ(-γ/λ) in virtue of the homeomorphism in [START_REF] Schmitz | Wasserstein dictionary learning: Optimal transport-based unsupervised non-linear dictionary learning[END_REF] and [START_REF] Schmitzer | A sparse multi-scale algorithm for dense optimal transport[END_REF].

Lemma 8. The restriction of the regularized cost •, γ + λφ(•) to the transport polytope Π(p, q)

attains its global minimum uniquely.

Proof. We notice that the regularized cost is equal to a Bregman divergence up to a positive factor and additive constant:

π, γ + λφ(π) -λφ(ξ) = λB φ (π ξ) . ( 64
)
Hence, its minimization over the closed convex set Π(p, q) is equivalent to the Bregman projection of ξ ∈ int(dom φ) onto Π(p, q) according to the function φ of Legendre type. Since Π(p, q) ∩ int(dom φ) = ∅, this projection exists and is unique.

Definition 2. The dual rot mover's distance is the quantity defined as:

d γ,λ,φ (p, q) = π λ , γ , ( 65
)
where the dual rot mover's plan π λ is given by:

π λ = argmin π∈Π(p,q) π, γ + λφ(π) . ( 66
)
Remark 5. For the sake of notation, we omit the dependence on p, q, γ, φ in the index of dual rot mover's plans π λ .

We proceed with the following proposition that enlightens the relation between the RMD and associated Bregman divergence.

Proposition 9. The dual rot mover's plan is the Bregman projection of ξ onto the transport polytope:

π λ = argmin π∈Π(p,q) B φ (π ξ) . ( 67
)
Proof. This is a consequence of the proof for Lemma 8. Indeed, from the definition in (66), we see that the rot mover's plan also minimizes (64). Therefore, it is the unique Bregman projection of ξ onto the transport polytope.

We have a geometrical interpretation where the regularization shrinks the solution toward the matrix ξ that has minimal Bregman information.

Proposition 10. The dual rot mover's plan π λ can be obtained as:

π λ = argmin π∈Π(p,q) π, γ + λB φ (π ξ ) . ( 68
)
Proof. Developing the Bregman divergence based on its definition [START_REF] Dhillon | Matrix nearness problems with Bregman divergences[END_REF], we have:

B φ (π ξ ) = φ(π) -φ(ξ ) -π -ξ , ∇φ(ξ ) . ( 69
)
Since ∇φ(ξ ) = 0, the last term with scalar product vanishes and we are left out with φ(π) plus a constant term with respect to π. Hence, we can replace φ(π) by B φ (π ξ ) in the minimization (66) that defines π λ .

Under some additional conditions, this interpretation can also be seen as shrinking toward the transport plan π with minimal Bregman information. Proposition 11. If π ∈ ri(Π(p, q)), then the dual rot mover's plan π λ can be obtained as:

π λ = argmin π∈Π(p,q) π, γ + λB φ (π π ) . ( 70
)
Proof. If π ∈ ri(Π(p, q)), then we have equality in the generalized Pythagorean theorem [START_REF] Haasdonk | Feature space interpretation of SVMs with indefinite kernels[END_REF], leading to:

B φ (π ξ ) = B φ (π π ) + B φ (π ξ ) . (71)
Since the last term is constant with respect to π, we can replace B φ (π ξ ) by B φ (π π ) in the minimization (68) that characterizes π λ .

Remark 6. The proposition also holds trivially when the global minimum is attained on the transport polytope, that is, when ξ = π .

Corollary 12. Under assumptions (A), the dual rot mover's plan π λ can be obtained as:

π λ = argmin π∈Π(p,q) π, γ + λB φ (π π ) . ( 72
)
Proof. This is a result of π ∈ Π(p, q) ∩ int(dom φ) = ri(Π(p, q)) when dom φ ⊆ R d×d + , as shown in the proof of Corollary 6.

In the sequel, we also extend naturally the definition of the dual RMD for λ = 0 as the EMD. We then do not necessarily have uniqueness of dual rot mover's plans for λ = 0, and the geometrical interpretation in terms of a Bregman projection does not hold anymore for λ = 0. However, we have the following theorem based on duality theory that shows the equivalence between primal and dual ROT problems.

Theorem 13. For all α > 0, there exists λ ≥ 0 such that the primal and dual rot mover's distances are equal:

d γ,α,φ (p, q) = d γ,λ,φ (p, q) . ( 73
)
Moreover, if α < α , then a corresponding value is such that λ > 0, and the primal and dual rot mover's plans are unique and equal:

π α = π λ . ( 74
)
Proof. The primal problem can be seen as the minimization p of the cost π, γ on Π(p, q) subject to φ(π) -φ(π ) -α ≤ 0. The domain of this constrained convex problem is D = Π(p, q)∩dom φ = ∅. The Lagrangian on D×R is given by L(π, λ) = π, γ +λ(φ(π)-φ(π )-α), and its minimization over D for a fixed λ ≥ 0 has the same solutions π as the dual problem.

In addition, Slater's condition for convex problems, stating that there is a strictly feasible point in the relative interior of the domain, is verified as long as α > 0. Indeed, we have ri(D) = ri(Π(p, q)). The existence of a strictly feasible point φ(π) < φ(π ) + α then holds by continuity of φ at π ∈ int(dom φ). As a result, we have strong duality with a zero duality gap p = d , where d is the maximization of g(λ) subject to λ ≥ 0. Moreover, if d is finite, then it is attained at least once at a point λ . This is the case since we already know that p is finite. Since p is also attained at least once at a point π solution of the primal problem, we have the following chain:

p = d (75) = min π∈D L(π, λ ) (76) ≤ L(π , λ ) (77) = π , γ + λ (φ(π ) -φ(π ) -α) (78) ≤ π , γ (79) = p . ( 80
)
Therefore, all inequalities are in fact equalities, π also minimizes the Lagrangian over D and thus is a solution of the dual problem. In other words, the primal and dual RMD for α and λ are equal, and the primal solutions must be dual solutions too. For α < α , the RMD has not reached the EMD yet, and thus we must have λ > 0. Hence, the dual solution is unique, so that the primal solution is unique too and equal to the dual one.

Remark 7. Corresponding values of α and λ depend on p, q, γ, φ. In addition, there might be multiple values of λ that correspond to a given α.

Again, the RMD does not verify the triangular inequality in general, and hence does not provide a true distance metric on Σ d even if γ is a distance matrix. Nevertheless, we still have the result that the RMD is symmetric as soon as φ is invariant by transposition, which holds for separable regularizers φ ij = φ, and γ is symmetric. We also obtain properties for the dual RMD that are similar to the ones for the primal RMD. d γ,λ,φ (p,q) is an increasing and continuous function of λ.

Property 6. The dual rot mover's distance

Proof. The fact that it is increasing can be seen as follows. Let 0 ≤ λ 1 < λ 2 . By construction, we have the following inequalities:

π λ1 , γ + λ 1 φ(π λ1 ) ≤ π λ2 , γ + λ 1 φ(π λ2 ) , ( 81
)
π λ2 , γ + λ 2 φ(π λ2 ) ≤ π λ1 , γ + λ 2 φ(π λ1 ) . ( 82
)
Subtracting these two inequalities, we obtain that φ(π λ1 ) ≥ φ(π λ2 ). Reinserting this result in the first inequality, we finally get π λ1 , γ ≤ π λ2 , γ . The continuity of the dual RMD results from that of the primal RMD. Let λ ≥ 0, and choose an arbitrary dual rot mover's plan π λ and earth mover's plan π . On the one hand, we have π , γ ≤ π λ , γ . On the other hand, we have

π λ , γ + λφ(π λ ) ≤ π , γ + λφ(π ), and thus π λ , γ ≤ π , γ + λ(φ(π ) -φ(π λ )) ≤ π , γ .
Suppose we have a discontinuity of the dual RMD at λ. Then by monotonicity, there is a value π , γ < d < π , γ that is not in the image of the dual RMD. But d is in the image of the primal RMD for a given α > 0 by continuity. It means that λ > 0 such that π λ , C = d, whereas, by continuity of the primal problem, we know that there exist α > 0 such that π α , C ≥ d. This is in contradiction with the duality result in Theorem 13, which implies that the image of the primal RMD for α > 0 must be included in that of the dual RMD for λ ≥ 0.

Property 7. When λ tends to +∞, the dual rot mover's distance converges to:

lim λ→+∞ d γ,λ,φ (p, q) = π , γ , ( 83
)
and the dual rot mover's plan converges in norm to the transport plan with minimal Bregman information:

lim λ→+∞ π λ = π . ( 84
)
Proof. Let (λ k ) k∈N be a sequence of positive numbers that tends to +∞, and (π λ k ) k∈N the associated rot mover's plans. By compactness of Π(p, q), we can extract a subsequence of rot mover's plans that converges in norm to a point π ∈ Π(p, q). For the sake of simplicity, we do not relabel this subsequence. By construction, we have

π λ k , γ + λ k φ(π ) ≤ π λ k , γ + λ k φ(π λ k ) ≤ π , γ +λ k φ(π ).
The scalar products are bounded, so dividing the inequalities by λ k and taking the limit, we obtain that φ(π λ k ) converges to φ(π ). By lower semi-continuity of φ, we thus have φ(π ) ≤ φ(π ). Since the global minimum of φ on Π(p, q) is attained uniquely at π , we must have π = π , and the original sequence also converges in norm to π . Hence, the dual rot mover's plan π λ converges in norm to π when λ tends to +∞. By continuity of the total cost •, γ on R d×d , π λ k , γ converges to π , γ . Hence, the limit of the RMD when λ tends to +∞ is π , γ . Property 8. When λ tends to 0, the dual rot mover's distance converges to the earth mover's distance:

lim λ→0 d γ,λ,φ (p, q) = d γ (p, q) . ( 85
)
Proof. This is a direct consequence of the dual RMD being continuous at λ = 0.

Property 9. If [0, 1) d×d ⊆ dom φ and φ is strictly convex on [0, 1) d×d , then the dual rot mover's plan converges in norm when λ tends to 0 to the earth mover's plan π 0 with minimal Bregman information:

lim λ→0 π λ = π 0 . ( 86
)
Proof. Let (λ k ) k∈N be a sequence of positive numbers that converges to 0, and (π λ k ) k∈N the associated rot mover's plans. By compactness of Π(p, q), we can extract a subsequence of rot mover's plans that converges in norm to a point π ∈ Π(p, q). For the sake of simplicity, we do not relabel this subsequence. By construction, we have

π 0 , γ + λ k φ(π λ k ) ≤ π λ k , γ + λ k φ(π λ k ) ≤ π 0 , γ + λ k φ(π 0 ).
The regularizer φ is continuous on the polytope Π(p, q) ⊆ dom φ, so taking the limit, we obtain that π λ k , γ converges to π 0 , γ . Therefore, π must be an earth mover's plan. Now dividing by λ k and taking the limit, we obtain that φ(π ) ≤ φ(π 0 ). Since π 0 is the unique earth mover's plan with minimal Bregman information, we must have π = π 0 .

Geometrical Insights

Our primal and dual formulations enlighten some intricate relations between optimal transportation theory [START_REF] Villani | Optimal Transport: Old and New[END_REF] and information geometry [START_REF] Amari | Methods of Information Geometry[END_REF], where Bregman divergences are known to possess a dually flat structure with a generalized Pythagorean theorem for information projections. A schematic view of the underlying geometry for ROT problems is represented in Figure 1, and can be discussed as follows.

Our constructions start from the global minimizer ξ of the regularizer φ (Lemma 1). The Bregman projection π of ξ onto the transport polytope Π(p, q) has minimal Bregman information on Π(p, q) (Lemma 2). The linear cost restricted to the regularized transport polytope Π α,φ (p, q) also attains its global minimum (Lemma 3). Such a minimizer π α is a primal rot mover's plan (Definition 1). We can interpret Π α,φ (p, q) as the intersection of Π(p, q) with the Bregman ball of radius B φ (π ξ ) + α and center ξ (Proposition 4). In certain cases, Π α,φ (p, q) is also the intersection of Π(p, q) with the Bregman ball of radius α and center π , as a result of the generalized Pythagorean theorem B φ (π ξ ) = B φ (π π ) + B φ (π ξ ) (Proposition 5, Corollary 6). All in all, this enforces the solutions to have small enough Bregman information, by constraining them to lie close to the matrix ξ or transport plan π with minimal Bregman information.

In our developments, we next introduce the global minimizer ξ of the regularized cost (Lemma 7). The regularized cost restricted to Π(p, q) also attains its global minimum uniquely (Lemma 8). This minimizer defines the dual rot mover's plan π λ (Definition 2). Actually, π λ can be seen as the Bregman projection of ξ onto Π(p, q) (Proposition 9). The regularization by the Bregman information is also equivalent to regularizing the solution toward ξ (Proposition 10). In some cases, this can also be seen as regularizing toward π , as a result of the generalized Pythagorean theorem B φ (π ξ ) = B φ (π π ) + B φ (π ξ ) (Proposition 11, Corollary 12). Again, Figure 1: Geometry of regularized optimal transport. this enforces the solutions to have small enough Bregman information, by shrinking them toward the matrix ξ or transport plan π with minimal Bregman information.

We have duality between the primal and dual formulations, so that primal and dual rot mover's plans follow the same path on Π(p, q) from no regularization (α = +∞, λ = 0) to full regularization (α = 0, λ = +∞) (Theorem 13). In the limit of no regularization, we obviously retrieve an earth mover's plan π for the cost matrix γ. By duality, it is also intuitive that the additional constraint for the primal formulation, seen in the equivalent forms of φ(π), B φ (π ξ ) or B φ (π π ), leads to an analog penalty for the dual formulation in the same respective form.

Since ξ = 1, π = pq , ξ = exp(-γ/λ), for minus the Boltzmann-Shannon entropy and Kullback-Leibler divergence, we retrieve the existing results discussed in Section 1.2 as a specific case [START_REF] Benamou | Iterative Bregman projections for regularized transportation problems[END_REF][START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF]. In addition, we can readily generalize the estimation of contingency tables with fixed marginals to a matrix nearness problem based on other divergences than the Kullback-Leibler divergence [START_REF] Dhillon | Matrix nearness problems with Bregman divergences[END_REF]. Given a rough estimate ξ ∈ int(dom φ), a contingency table with fixed marginals p, q can be estimated by Bregman projection of ξ onto Π(p, q):

π = argmin π∈Π(p,q) B φ (π ξ) . ( 87
)
This simply amounts to solving a dual ROT problem with an arbitrary penalty λ > 0 and a cost matrix γ = -λ∇φ(ξ).

Finally, since Bregman divergences are invariant under adding an affine term to their generator, it is straightforward to generalize ROT problems by shrinking toward an arbitrary prior matrix ξ, ξ ∈ int(dom φ), or transport plan π ∈ Π(p, q). This is indeed equivalent to translating the regularizer by the appropriate amount φ(π) + π, δ , so that the global minimizer is now attained at the desired point. Equivalently, this amounts to translating the cost matrix as γ + λδ instead.

Algorithmic Derivations

In this section, we introduce algorithmic methods to solve ROT problems. We focus without lack of generality on the dual problem, which can be solved efficiently via alternate Bregman projections. The primal problem can then easily be solved for 0 < α < α by a bisection search on λ > 0. For α = 0, we could simply use alternate Bregman projections to project ∇ψ(0) instead of ∇ψ(-γ/λ) in virtue of Lemmas 1 and 2, which actually corresponds to the special case γ = 0 in our algorithms, though this is not really relevant in practice since this completely removes the linear influence of the total cost from the ROT problem. In the limit α ≥ α , a classical OT solver such as the network simplex can directly be used. We first study the underlying Bregman projections in their generic form (Section 4.1) and specifically develop the case of separable divergences (Section 4.2). We then derive the two generic schemes of ASA (Section 4.3) and NASA (Section 4.4) to solve dual ROT problems, depending on whether the domain of the smooth convex regularizer lies within the non-negative orthant or not. We also enhance both algorithms in the separable case with a sparse extension (Section 4.5), and finally discuss some practical considerations of our methods (Section 4.6). To simplify notations, we omit the penalty value λ in the index and simply write π for the rot mover's plan.

Generic Projections

The closed convex transport polytope Π(p, q) is the intersection of the non-negative orthant:

C 0 = R d×d + , (88)
which is a polyhedral subset, with two affine subspaces:

C 1 = {π ∈ R d×d : π1 = p} , ( 89
)
C 2 = {π ∈ R d×d : π 1 = q} . ( 90
)
The Bregman projection π onto Π(p, q) can then be obtained by alternate Bregman projections onto C 0 , C 1 , C 2 , where we expect that these latter projections are easier to compute. On the one hand, the Karush-Kuhn-Tucker conditions for Bregman projection π 0 of a given matrix π ∈ int(dom φ) onto C 0 are necessary and sufficient, and write as follows:

π 0 ≥ 0 , ( 91
)
∇φ(π 0 ) -∇φ(π) ≥ 0 , ( 92
) (∇φ(π 0 ) -∇φ(π)) π 0 = 0 . ( 93
)
While these conditions are nontrivial to solve in general, we shall see that they admit an elegant solver specific to the non-separable squared Mahalanobis distances defined in [START_REF] Montavon | Wasserstein training of restricted Boltzmann machines[END_REF] and generated by the quadratic form in [START_REF] Muzellec | Tsallis regularized optimal transport and ecological inference[END_REF]. In addition, they also greatly simplify for separable divergences, which encompass all other divergences used in this paper.

On the other hand, the Lagrangians with Lagrange multipliers µ, ν ∈ R d for the Bregman projections π 1 and π 2 of a given matrix π ∈ int(dom φ) onto C 1 and C 2 respectively write as follows:

L 1 (π, µ) = φ(π) -π, ∇φ(π) + µ (π1 -p) , ( 94
) L 2 (π, ν) = φ(π) -π, ∇φ(π) + ν (π 1 -q) . ( 95
)
Their gradients are given on int(dom φ) by:

∇L 1 (π, µ) = ∇φ(π) -∇φ(π) + µ1 , ( 96
)
∇L 2 (π, ν) = ∇φ(π) -∇φ(π) + 1ν , ( 97
)
and vanish at π 1 , π 2 ∈ int(dom φ) if and only if:

π 1 = ∇ψ(∇φ(π) -µ1 ) , ( 98
)
π 2 = ∇ψ(∇φ(π) -1ν ) . ( 99
)
By duality, the Bregman projections onto C 1 , C 2 are thus equivalent to finding the unique vectors µ, ν, such that the rows of π 1 sum up to p, respectively the columns of π 2 sum up to q:

∇ψ(∇φ(π) -µ1 )1 = p , ( 100
)
∇ψ(∇φ(π) -1ν ) 1 = q . ( 101
)
Similarly, solving for the Lagrange multipliers is an expensive problem in general, since the search space is of dimension d and we evaluate matrix functions of size d × d. This is because a given entry µ i , ν j can actually modify any entry of the d × d matrix functions being evaluated. Again, we shall see that they can nevertheless be computed efficiently for separable divergences as well as the non-separable Mahalanobis distances.

Separable Case

Assuming that the regularizer φ is separable, the underlying Bregman projections can be computed more efficiently. To keep notations simple, we focus on separable divergences with same element-wise regularizer, and thus chiefly omit the indices φ ij = φ. We emphasize, however, that it is straightforward to apply all our methods for separable divergences with different elementwise regularizers, which notably enables weighting a given element-wise regularizer.

In case of separability, the Karush-Kuhn-Tucker conditions for projection onto C 0 simplify to provide a closed-form solution on primal parameters:

π 0,ij = max{0, π ij } . ( 102
)
Since φ is increasing, this is equivalent on dual parameters to:

θ 0,ij = max{φ (0), θ ij } . ( 103
)
Now turning to projections onto C 1 , C 2 for primal parameters π 1,ij , π 2,ij , we can divide the initial problems into d parallel subproblems in search space of dimension 1 each. This is much more efficient to solve than in the non-separable case. This can be summarized as looking for d separate Lagrange multipliers µ i , respectively ν j , such that:

d j=1 ψ (θ ij -µ i ) = p i , ( 104
) d i=1 ψ (θ ij -ν j ) = q j . ( 105
)
Finding the optimal values µ i , ν j ∈ R through ψ and the sums over rows or columns, however, is still nontrivial in general.

An analytical solution can be obtained in specific cases. Intuitively, we need to factor µ i , ν j out of ψ as additive or multiplicative terms. This is related to Pexider's functional equations, which hold only for functions with a linear form ψ (θ) = aθ + b, or exponential form ψ (θ) = a exp(bθ), with a, b ∈ R. This leads to regularizers with a quadratic form φ(π) = aπ 2 + bπ + c, or entropic form φ(π) = aπ log π +bπ +c, with a, b, c ∈ R. The constants a, b actually only scale and translate the cost matrix, whereas the constant c has no effect. Referring to Table 1, the quadratic case holds under assumptions (B), and thus requires Dykstra's algorithm for alternate Bregman projections with correction terms to ensure non-negativity by projection onto the polyhedral nonnegative orthant. The entropic case holds under assumptions (A), using the POCS technique for alternate Bregman projection with no correction terms since the non-negativity is already ensured by the domain of the regularizer. The latter case reduces to the regularization of [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] and [START_REF] Benamou | Iterative Bregman projections for regularized transportation problems[END_REF], so that we actually end up with the Sinkhorn-Knopp algorithm. Hence, the Euclidean norm associated to the squared Euclidean distance, and the entropic case associated to the Kullback-Leibler divergence, are reasonably the only two existing analytical schemes to find the sum constraint projections. For other ROT problems, available solvers for line search can be employed instead.

For simplicity, we assume hereafter that ψ is twice continuously differentiable with ψ positive and ψ verifying the necessary and sufficient condition [START_REF] Rolet | Fast dictionary learning with a smoothed Wasserstein loss[END_REF] on its whole domain. Therefore, we can use the Newton-Raphson method with guarantees of global convergence. This encompasses most of the common regularizers, and notably all regularizers used in this paper except from the Fermi-Dirac entropy, p norms and Hellinger distance. When the condition (40) for global convergence is not met on the whole domain, it is still possible to apply the Newton-Raphson method after careful initialization, so as to restrict to a smaller interval where the condition holds. This is discussed in more detail with practical examples for the Fermi-Dirac entropy, p norms and Hellinger distance in Section 5, where the first-order derivatives are increasing convex on half of the domain and increasing concave on the other half. When the second-order derivatives do not exist, are not continuous or vanish at some points, a similar strategy can be applied. This is again discussed for the p norms in Section 5, where the second-order derivative is undefined or vanishes at 0 depending on the value of the parameter. If such an initialization is not possible, then a bisection search can always be applied instead of the Newton-Raphson method.

To apply the Newton-Raphson method, we exploit the following functions:

f (µ i ) = - d j=1 ψ (θ ij -µ i ) , ( 106
) g(ν j ) = - d i=1 ψ (θ ij -ν j ) , ( 107
)
defined respectively on the open intervals ( θi -θ, +∞) and ( θj -θ, +∞), where 0 < θ ≤ +∞ is such that dom ψ = (-∞, θ), and θi = max{θ ij } 1≤j≤d , θj = max{θ ij } 1≤i≤d . Their continuous derivatives are given by:

f (µ i ) = d j=1 ψ (θ ij -µ i ) , ( 108
) g (ν j ) = d i=1 ψ (θ ij -ν j ) , (109)
and are positive, so that f, g are strictly increasing on their whole domain, and thus on any closed interval with endpoints consisting of a feasible point and a solution. By construction, f, g also verify the necessary and sufficient condition [START_REF] Rolet | Fast dictionary learning with a smoothed Wasserstein loss[END_REF] for global convergence, and we know that there are unique solutions to f (µ i ) = -p i and g(ν j ) = -q j . Hence, the Newton-Raphson updates:

µ i ← µ i + d j=1 ψ (θ ij -µ i ) -p i d j=1 ψ (θ ij -µ i ) , ( 110
)
ν j ← ν j + d i=1 ψ (θ ij -ν j ) -q j d i=1 ψ (θ ij -ν j ) , ( 111
)
converge to the optimal solutions with a quadratic rate for any feasible starting points. By construction, we also know that initialization can be done with µ i ← 0, ν j ← 0. To avoid storing the intermediate Lagrange multipliers, the updates can then directly be written on dual parameters:

θ 1,ij ← θ 1,ij - d j=1 ψ (θ 1,ij ) -p i d j=1 ψ (θ 1,ij ) , ( 112
)
θ 2,ij ← θ 2,ij - d i=1 ψ (θ 2,ij ) -q j d i=1 ψ (θ 2,ij ) , ( 113
)
after initialization by θ 1,ij ← θ ij , θ 2,ij ← θ ij .

Alternate Scaling Algorithm

Under assumptions (A), we can drop the non-negative constraint since it is already ensured by dom φ ⊆ R d×d + (Table 1). The POCS technique in its basic form [START_REF] Pele | Fast and robust earth mover's distances[END_REF] then states that the projection of ξ onto Π(p, q) can be obtained by alternate Bregman projections onto the affine subspaces C 1 and C 2 with linear convergence. Clearly, the underlying control mapping takes each output value an infinite number of times. Since we have just two sets, the only possible alternative in the control mapping is to swap the order of projections starting from C 2 instead of C 1 , which actually amounts to swapping the input distributions p, q and transposing the cost matrix γ, to obtain the transposed of the rot mover's plan. We thus focus on the first choice without lack of generality.

Starting from ξ and writing the successive vectors µ (k) , ν (k) along iterations, we have the following sequence:

∇ψ(-γ/λ) → ∇ψ -γ/λ -µ (1) 1 (114)
→ ∇ψ -γ/λ -µ (1) 1 -1ν (1) (115)

→ . . . ( 116
)
→ ∇ψ -γ/λ -µ (1) 1 -1ν (1)

-• • • -µ (k) 1 (117) → ∇ψ -γ/λ -µ (1) 1 -1ν (1) -• • • -µ (k) 1 -1ν (k) (118)
Algorithm 1 Alternate scaling algorithm. θ ← -γ/λ repeat θ ← θ -µ1 , where µ uniquely solves ∇ψ(θ -µ1 )1 = p θ ← θ -1ν , where ν uniquely solves ∇ψ(θ -1ν ) 1 = q until convergence π ← ∇ψ(θ ) Algorithm 2 Alternate scaling algorithm in the separable case.

θ ← -γ/λ repeat repeat θ ← θ -ψ (θ )1-p ψ (θ )1 1 until convergence repeat θ ← θ -1 1 ψ (θ )-q 1 ψ (θ ) until convergence until convergence π ← ψ (θ ) → . . . ( 119
) → π . ( 120
)
In other terms, we obtain the rot mover's plan π by scaling iteratively the rows and columns of the successive estimates through ∇ψ. An efficient algorithm, called ASA, is to store a unique d × d matrix in dual parameter space and update it by alternating the projections in primal parameter space (Algorithm 1). The updates have a complexity in O(d 2 ) once the vectors µ, ν are obtained.

In the separable case, the projections can be obtained by iterating the respective Newton-Raphson update steps, which can be written compactly with matrix and vector operations (Algorithm 2). The complexity for the updates are now clearly in O(d 2 ). In more detail, each update step features one vector row or column replication, one vector element-wise division, one vector subtraction, one matrix subtraction, two matrix row or column sums, and two elementwise matrix function evaluations. Because of separability, we can expect the required number of iterations for convergence in the different loops to be independent of the data dimension, and thus expect a quadratic empirical complexity as well.

Non-negative Alternate Scaling Algorithm

Under assumptions (B), we must now include the non-negative constraint since dom φ R d×d + (Table 1). We suggest to ensure non-negativity of each update, and thus follow a cycle of projections onto C 0 , C 1 , C 0 , C 2 . The underlying control mapping is a fortiori essentially cyclic. For practical reasons, we also ensure non-negativity of the output solution with a final projection onto C 0 . Again, swapping the order of projections onto C 1 , C 2 is equivalent to swapping the input distributions p, q and transposing the cost matrix γ to obtain the transposed of the rot mover's plan. Other control mappings could also be exploited, for example by ensuring non-negativity every two or more sum constraint projections. We do not discuss such variants here and focus Algorithm 3 Non-negative alternate scaling algorithm.

θ ← -γ/λ ϑ ← 0 ← 0 ς ← 0 θ ← θ + ϑ θ ← θ, where θ uniquely solves ∇ψ(θ) ≥ 0, θ ≥ θ, (θ -θ) ∇ψ(θ) = 0 ϑ ← θ -θ repeat θ ← θ + θ ← θ -µ1 , where µ uniquely solves ∇ψ(θ -µ1 )1 = p ← θ -θ θ ← θ + ϑ θ ← θ, where θ uniquely solves ∇ψ(θ) ≥ 0, θ ≥ θ, (θ -θ) ∇ψ(θ) = 0 ϑ ← θ -θ θ ← θ + ς θ ← θ -1ν , where ν uniquely solves ∇ψ(θ -1ν ) 1 = q ς ← θ -θ θ ← θ + ϑ θ ← θ, where θ uniquely solves ∇ψ(θ) ≥ 0, θ ≥ θ, (θ -θ) ∇ψ(θ) = 0 ϑ ← θ -θ until convergence π ← ∇ψ(θ )
on the above-mentioned sequence. The non-negative orthant being polyhedral but not affine, we also need to incorporate correction terms ϑ, , ς for all three projections. In more detail, the projections are computed after correction so that we do not directly project the obtained updates θ but the corrected updates θ = θ + ϑ, θ = θ + , and θ = θ + ς for the respective subsets. The correction terms are also updated as the difference θθ between the projected point and its projection. Dykstra's algorithm [START_REF] Oberman | An efficient linear programming method for optimal transportation[END_REF] for Bregman divergences with corrections [START_REF] Pele | A linear time histogram metric for improved SIFT matching[END_REF] then guarantees that the projection of ξ onto Π(p, q) is obtained with linear convergence.

A general algorithm, called NASA, is to store d × d matrices for projected points, projections and correction terms in dual parameter space, update them accordingly and finally go back to primal parameter space (Algorithm 3). The updates have a complexity in O(d 2 ) once the Karush-Kuhn-Tucker conditions are solved or Lagrange multipliers µ, ν are obtained.

In the separable case, the non-negativity constraint can be obtained analytically and the sequence of updates greatly simplifies. Starting from ξ and writing the successive vectors µ (k) , ν (k) along iterations, we have:

ψ (-γ/λ) → ψ max{φ (0), -γ/λ}
→ ψ max{φ (0), -γ/λ} -µ (1) 1 → ψ max{φ (0), -γ/λ -µ (1) 1 } → ψ max{φ (0), -γ/λ -µ (1) 1 } -1ν (1) → ψ max{φ (0), -γ/λ -µ (1) 1 -1ν (1) } Algorithm 4 Non-negative alternate scaling algorithm in the separable case.

θ ← -γ/λ θ ← max{φ (0), θ} repeat τ ← 0 repeat τ ← τ + ψ (θ -τ 1 )1-p ψ (θ -τ 1 )1 until convergence θ ← θ -τ 1 θ ← max{φ (0), θ} σ ← 0 repeat σ ← σ + 1 ψ (θ -1σ )-q 1 ψ (θ -1σ ) until convergence θ ← θ -1σ θ ← max{φ (0), θ} until convergence π ← ψ (θ )
→ ψ max{φ (0), -γ/λ -µ (1) 1 -1ν (1) } + µ (1) 1 -µ (2) 1 → ψ max{φ (0), -γ/λ -µ (2) 1 -1ν (1) } → ψ max{φ (0), -γ/λ -µ (2) 1 -1ν (1) } + 1ν (1) -1ν (2) → ψ max{φ (0), -γ/λ -µ (2) 1 -1ν (2) } → . . .

→ ψ max{φ (0), -γ/λ -µ (k) 1 -1ν (k) } → ψ max{φ (0), -γ/λ -µ (k) 1 -1ν (k) } + µ (k) 1 -µ (k+1) 1 → ψ max{φ (0), -γ/λ -µ (k+1) 1 -1ν (k) } → ψ max{φ (0), -γ/λ -µ (k+1) 1 -1ν (k) } + 1ν (k) -1ν (k+1) → ψ max{φ (0), -γ/λ -µ (k+1) 1 -1ν (k+1) } → . . . → π .
An efficient algorithm then exploits the differences τ (k) = µ (k) -µ (k-1) and σ (k) = ν (k) -ν (k-1) to scale the rows and columns (Algorithm 4). We store d×d matrices as well as difference vectors instead of correction matrices. The algorithm can then be interpreted as producing interleaved updates between the projections according to the max operator and according to the respective scalings. The updates in NASA now clearly have a complexity in O(d 2 ) when using the Newton-Raphson method for scaling, with similar matrix and vector operations to ASA in the separable case, and an expected empirical complexity that is quadratic.

Sparse Extension

In the separable case, it is possible to develop a sparse extension of both our methods ASA and NASA. Storing and updating full d × d matrices becomes expensive with the data dimension. Instead, we allow for infinite entries in the cost matrix γ, meaning that the transport of mass between certain bins is proscribed. As a result, the corresponding entries of π must be null. Eventually, we can drop all these entries so that we just need to store and update the remaining ones. The RMD via the Frobenius inner product π , γ is then computed without accounting for discarded entries, or equivalently by setting indefinite element-wise products 0 × ∞ = 0 by convention, so it naturally costs nothing to move no mass on a path that is forbidden. This leads to an expected complexity in O(r), where r is the number of finite entries in γ. Typically, r can be chosen in the order of magnitude of d, so as to obtain a linear instead of quadratic empirical complexity.

In practice, both ASA and NASA are compatible with this strategy. We always have lim θ→-∞ ψ (θ) = 0 under assumptions (A) for ASA. Under assumptions (B) for NASA, this limit might be finite or infinite but is necessarily negative, so also leads to 0 after enforcing non-negativity by projection onto the non-negative orthant. As a result, the obtained sequence of projections preserves the desired zeros in both algorithms, and an infinite element-wise cost does lead to no mass transport at all between the corresponding bins. In theory, we can understand this extension in light of the dual formulation seen as a Bregman projection in (67). Under assumptions (B), we always have 0 ∈ int(dom φ) and thus B φ (0 0) = 0. Hence, Dykstra's algorithm is readily applicable in the sparse version. Under assumptions (A), however, we have 0 / ∈ int(dom φ), and even sometimes 0 / ∈ dom φ as for the Itakura-Saito divergence. We can nonetheless extend the domain of the element-wise divergence at the origin by continuity on the diagonal, that is, by setting it null as B φ (0 0) = 0. This is akin to considering absolutely continuous measures, also known as dominated measures, and Radon-Nikodym derivatives to generalize the definition of Bregman divergences. [START_REF] Kurras | Symmetric iterative proportional fitting[END_REF] then showed that the POCS method still holds with this convention by introducing a notion of locally affine spaces.

With such a sparse extension, however, we must take care that a sparse solution does exist, meaning that there is a transport plan in the transport polytope that has the desired zeros. For example, if all entries of γ are infinite, then there are obviously no possible sparse solutions since we enforce all entries of the plan to be null. A necessary condition for the existence of a sparse solution is that for any entry q j , all entries p k from which we are allowed to transport mass must provide enough total mass to fill q j completely. Similarly, for any entry p i , all entries q k to which we are allowed to transport mass must require enough total mass to empty p i completely. Unfortunately, sufficient conditions are not so intuitive. [START_REF] Idel | A review of matrix scaling and Sinkhorn's normal form for matrices and positive maps[END_REF]Theorem 4.1] studied such problems thoroughly and elucidated several necessary and sufficient conditions for sparse solutions to exist, but these conditions are nontrivial to use from in practice. [START_REF] Kurras | Symmetric iterative proportional fitting[END_REF] advocates trying first to compute a solution with the desired sparsity, and if no solution can be found, then gradually reduce sparsity until a solution is found. This might still speed up computation drastically because of the linear instead of quadratic complexity. Lastly, we remark that it is not evident to propose a sparse extension for the non-separable case in general, since a given entry of γ might influence all entries of π .

Practical Considerations

As noticed by [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] and [START_REF] Benamou | Iterative Bregman projections for regularized transportation problems[END_REF], the Sinkhorn-Knopp algorithm might fail to converge because of numerical instability when the penalty λ gets small. In particular, unless taking special care of numerical stabilization [START_REF] Schmitzer | Stabilized sparse scaling algorithms for entropy regularized transport problems[END_REF], a direct limitation is the machine precision under which some entries of exp(-γ/λ) are represented as zeros in memory. Such issues occur similarly for other regu-larizations, notably via the representation ∇ψ(-γ/λ) of the unconstrained solution to project. Therefore, the proposed methods are actually competitive in a range where the penalty λ is not too small, and for which the rot mover's plan π exhibits a significant amount of smoothing. Hence, we do not target the same problems as traditional schemes such as interior point methods or the network simplex.

In addition, the different Bregman projections in our algorithms are most of the time approximate up to a given tolerance depending on the termination criterion used for convergence. Exceptions occur for the sum constraints with the Euclidean distance or Kullback-Leibler divergence, as well as the non-negativity constraints in the separable case, which are obtained analytically. A natural question to raise is then whether our algorithms still converge when the projections are approximate only. However, this is relatively hard to answer in theory. We did not observe in practice any problem of convergence when using sufficiently good approximations. Furthermore, first approximations can be quite rough without affecting convergence as long as final approximations are good enough. Sometimes, even alternating a single or two steps of the Newton-Raphson method throughout the main iterations the algorithm still works, though this is not systematic. Thus, we advocate for safety to use a tight tolerance for the auxiliary projections.

We also observed numerical instability of the Newton-Raphson updates for separable divergences under assumptions (A). This is due to the denominator being based on ψ with limit lim θ→-∞ ψ (θ) = 0, that is, for entries π close to zero. It is possible, however, to make the updates of µ i , ν j much more stable in practice by using the max truncation operator, despite theoretical guarantees of convergence without it. Specifically, we know that the entries π 1,ij must lie between 0 and p i , and π 2,ij between 0 and q j . Hence, we can lower bound µ i and ν j by θi -φ (p i ) and θj -φ (q j ), respectively. Interestingly, this also speeds up the convergence of the updates significantly when the initialization by 0 is far from the actual solution.

A possible termination criterion for the main and auxiliary iterations is to compute the marginal difference between the updated matrix and p, q. In the auxiliary iterations for the two scaling projections, we compare the sums of rows or columns to p or q respectively, and in the main iterations of the algorithm, we compare both marginals simultaneously. Typically, we can use the p (quasi-)norm with 0 < p ≤ +∞ to assess the marginal difference, and the auxiliary tolerance should be at least the main one for sufficient precision of the approximations. Two alternative quantities in absolute or relative scales can also be used for termination, either the variation with p (quasi-)norm in the updated matrix or in the updated distance. Here the auxiliary tolerance should be at least the square of the main one. This seems reasonable for π given the quadratic rate of convergence for the Newton-Raphson method versus the linear one for alternate Bregman projections, as well as for π , γ under the Cauchy-Schwarz inequality. In all cases, convergence can be checked either after each iteration or after a given number of iterations to reduce the underlying cost of computing the termination criterion. We can also fix a maximum number of main and auxiliary iterations to limit the overall running time.

Regarding implementation, the matrix and vector operations used for ASA and NASA in the separable case are well-suited for fast calculation on a GPU and for processing of multiple input distributions in parallel. By working directly in the primal parameter space, the Sinkhorn-Knopp algorithm is also readily suited for dealing with sparse plans, based on existing libraries. In more general ROT problems, however, a specific library should be written for the sparse extension because null entries in the transport plan are not represented by null entries in the dual parameter space, so that tailored data structures and operations for such matrices need to be coded. Therefore, we only implemented and will focus in our experiments on the non-sparse version of our methods.

Finally, although we implicitly assumed throughout that the entries of p and q are strictly -

β = 0 β = 0.5 β = 1 p = 0.1 p = 0.5 p = 0.9 p = 1.1 p = 1.5 p = 2 - FDLOG BETA BETA BETA LPQN LPQN LPQN LPN LPN LPN HELL - BIS - BSKL - - - - - EUC -
Figure 2: Separable regularizers on (0, 1).

comprised between 0 and 1 for theoretical issues, it is often possible in practice to deal explicitly with null or unit entries in the input distributions. Intuitively, no mass can be moved from or to a null entry, so the transport plans have null rows and columns for the corresponding null entries of p and q, respectively. In the separable case, we can thus simply remove these entries, solve the reduced ROT problem, and reinsert the corresponding null entries in the rot mover's plan π . The same reasoning as for the sparse extension can be made to show that our two algorithms still hold with this strategy from a theoretical standpoint. In the non-separable case, however, this is not as straightforward again because the influences of the different entries of π are interleaved through the regularizer φ. Nonetheless, as long as we have [0, 1) d×d ⊂ int(dom φ) under assumptions (B), then we have Π(p, q) ⊂ int(dom φ) and we can apply NASA without modification. This is notably the case for the Mahalanobis distances whose domain is R d×d . For a non-separable regularizer under assumptions (A), it is not easy to account for null entries because the constraint qualification Π(p, q) ∩ int(dom φ) = ∅ never holds due to mandatory null entries in the transport plans. Nevertheless, common regularizers under assumptions (A), including the ones used in this paper, are separable in general. Lastly, it is direct to cope with unit entries in p or q in all cases, since the transport polytope then reduces to a singleton, so that there is a unique transport plan pq which is the rot mover's plan.

Classical Regularizers and Divergences

In this section, we discuss the specificities of the ASA (Algorithm 1) and NASA (Algorithm 3) methods to solve ROT problems for classical regularizers and associated divergences. We start with several separable regularizers under assumptions (A), based on the Boltzmann-Shannon entropy related to the Kullback-Leibler divergence (BSKL, Section 5.1), the Burg entropy related to the Itakura-Saito divergence (BIS, Section 5.2), and the Fermi-Dirac entropy related to a logistic loss function (FDLOG, Section 5.3), as well as the parametric families of β-potentials related to the β-divergences (BETA, Section 5.4). We then discuss the separable p quasi-norms (LPQN, Section 5.5), which require a slight adaptation of assumptions (A). We also consider separable regularizers under assumptions (B) related to p norms (LPN, Section 5.6), as well as the Euclidean norm related to the Euclidean distance (EUC, Section 5.7) and the Hellinger distance (HELL, Section 5.8). Finally, we study a non-separable regularizer under assumptions (B) via quadratic forms in relation to Mahalanobis distances (Section 5.9). We plot all separable regularizers in Figure 2. All regularizers and their corresponding divergences are also summed up in Table 2. Lastly, we provide in Table 3 the related terms based on derivatives that are needed to instantiate the separable versions of ASA (Algorithm 2) or NASA (Algorithm 4) accordingly.

φ(π) / φ(π) B φ (π ξ) / B φ (π ξ) dom φ dom ψ Boltzmann-Shannon entropy Kullback-Leibler divergence π log π -π + 1 π log π ξ -π + ξ R + R Burg entropy Itakura-Saito divergence π -log π -1 π ξ -log π ξ -1 R ++ (-∞, 1)
Fermi-Dirac entropy Logistic loss function

π log π + (1 -π) log(1 -π) π log π ξ + (1 -π) log 1-π 1-ξ [0, 1] R β-potentials (0 < β < 1) β-divergences 1 β(β-1) (π β -βπ + β -1) 1 β(β-1) (π β + (β -1)ξ β -βπξ β-1 ) R + (-∞, 1 1-β ) p quasi-norms (0 < p < 1) -π p -π p + pπξ p-1 -(p -1)ξ p R + R -- p norms (1 < p < +∞) |π| p |π| p -pπ sgn(ξ)|ξ| p-1 + (p -1)|ξ| p R R Euclidean norm Euclidean distance 1 2 π 2 1 2 (π -ξ) 2 R R Hellinger distance -(1 -π 2 ) 1 2 (1 -πξ)(1 -ξ 2 ) -1 2 -(1 -π 2 ) 1 2 [-1, 1] R Quadratic forms (P 0) Mahalanobis distances 1 2 vec(π) P vec(π) 1 2 vec(π -ξ) P vec(π -ξ) R d×d R d×d
Table 2: Convex regularizers and associated Bregman divergences.

Boltzmann-Shannon Entropy and Kullback-Leibler Divergence

Assumptions (A) hold for minus the Boltzmann-Shannon entropy π log π -π + 1 associated to the Kullback-Leibler divergence. Hence, the ROT problem can be solved with the ASA scheme.

In addition, the updates in the POCS technique can be written analytically, leading to the Sinkhorn-Knopp algorithm. Specifically, the two projections amount to normalizing in turn the rows and columns of π so that they sum up to p and q respectively:

π ← diag p π 1 π , ( 121
) π ← π diag q π 1 . ( 122
)
This can be optimized by remarking that the iterates π (k) after each couple of projections verify:

π (k) = diag(u (k) )ξ diag(v (k) ) , ( 123
)
where ξ = exp(-γ/λ), and vectors u (k) , v (k) satisfy the following recursion:

u (k) = p ξv (k-1) , ( 124
)
v (k) = q ξ u (k) , ( 125
) φ(π) φ (π) ψ (θ) ψ (θ) Boltzmann-Shannon entropy π log π -π + 1 log π exp θ exp θ Burg entropy π -log π -1 1 -π -1 (1 -θ) -1 (1 -θ) -2
Fermi-Dirac entropy

π log π + (1 -π) log(1 -π) log π 1-π exp θ (1+exp θ) exp θ (1+exp θ) 2 β-potentials (0 < β < 1) 1 β(β-1) (π β -βπ + β -1) 1 β-1 (π β-1 -1) ((β -1)θ + 1) 1 β-1 ((β -1)θ + 1) 1 β-1 -1 p quasi-norms (0 < p < 1) -π p -pπ p-1 p -1 p-1 (-θ) 1 p-1 -p -1 p-1 p-1 (-θ) 1 p-1 -1 p norms (1 < p < +∞) |π| p p sgn(π)|π| p-1 p -1 p-1 sgn(θ)|θ| 1 p-1 p -1 p-1 p-1 |θ| 1 p-1 -1
Euclidean norm

1 2 π 2 π θ 1 Hellinger distance -(1 -π 2 ) 1 2 π(1 -π 2 ) -1 2 θ(1 + θ 2 ) -1 2 (1 + θ 2 ) -3 2
Table 3: Separable regularizers and related terms based on derivatives.

with convention v (0) = 1. This allows a fast implementation by performing only matrix-vector multiplications using a fixed matrix ξ = exp(-γ/λ). We can further save one element-wise vector multiplication per update:

u ← 1 diag 1 p ξ v , ( 126
) v ← 1 diag 1 q ξ u , ( 127
)
where the matrices diag 1 p ξ and diag 1 q ξ are precomputed and stored.

Burg Entropy and Itakura-Saito Divergence

Assumptions (A) also hold for minus the Burg entropy π -log π -1 associated to the Itakura-Saito divergence, so the ROT problem can be solved with the ASA scheme. Eventually, the Newton-Raphson steps to update the alternate projections in POCS technique can be written as follows:

θ ← θ - (1 -θ ) -1 1 -p (1 -θ ) -2 1 1 , ( 128
) θ ← θ -1 1 (1 -θ ) -1 -q 1 (1 -θ ) -2 . ( 129
)
Each step can be optimized by computing first an element-wise matrix inverse (1 -θ ) -1 for the numerator, and then performing an element-wise matrix multiplication of this matrix by itself to obtain a matrix for the denominator instead of applying an additional element-wise matrix power. Since ψ is convex and strictly increasing with ψ positive everywhere, the convergence of the updates is guaranteed.

Fermi-Dirac Entropy and Logistic Loss Function

Assumptions (A) again hold for minus the Fermi-Dirac entropy π log π + (1 -π) log(1 -π), also known as bit entropy, associated to a logistic loss function. The ROT problem can thus be solved with the ASA scheme, and the Newton-Raphson steps to update the alternate projections in the POCS technique can be written as follows:

θ ← θ - exp θ 1+exp θ 1 -p exp θ (1+exp θ ) 2 1 1 , ( 130
) θ ← θ -1 exp θ 1+exp θ -q 1 exp θ (1+exp θ ) 2 . ( 131
)
Each step can be optimized by storing first the element-wise matrix exponential exp θ , then applying an element-wise matrix division by the temporary matrix 1 + exp θ to obtain a matrix for the numerator, and lastly performing an element-wise matrix division of these two matrices to obtain a matrix for the denominator and thus save an additional element-wise matrix power as well as several element-wise matrix exponentials. However, even if ψ is strictly increasing with ψ positive everywhere, ψ is neither convex nor concave and does not verify the necessary and sufficient condition [START_REF] Rolet | Fast dictionary learning with a smoothed Wasserstein loss[END_REF] for global convergence of the Newton-Raphson method.

Nevertheless, ψ is convex on R -and concave on R + . It thus divides for a given 1

≤ i ≤ d, respectively 1 ≤ j ≤ d, the real line into at most d + 1 intervals -∞ < θ(1) i ≤ θ(2) i ≤ • • • ≤ θ(d-1) i ≤ θ(d) i < +∞, respectively -∞ < θ(1) j ≤ θ(2) j ≤ • • • ≤ θ(d-1) j ≤ θ(d) j < +∞, with
the values ( θ(k) i ) 1≤k≤d from row i of θ, respectively ( θ(k) j ) 1≤k≤d from column j of θ, sorted in increasing order. On each of these intervals, the necessary and sufficient condition ( 40) is verified since we can decompose f (µ i ), respectively g(ν j ), as the sum of an increasing convex and an increasing concave function. Hence, we have global convergence on the interval that contains the solution. It is further possible to restrict the search to the two last intervals only. Indeed, we have

d j=1 ψ (θ ij - θ(d-1) i ) ≥ ψ ( θ(d-1) i - θ(d-1) i ) + ψ ( θ(d) i - θ(d-1) i ) ≥ 2ψ (0) = 1, so that θ(d-1) i < µ i < +∞. Similarly, we have d i=1 ψ (θ ij - θ(d-1) j ) ≥ ψ ( θ(d-1) j - θ(d-1) j ) + ψ ( θ(d) j - θ(d-1) j ) ≥ 2ψ (0) = 1, so that θ(d-1) j < ν j < +∞. As a result, it suffices to initialize µ i with θi = θ(d) i = max {θ ij } 1≤j≤d , respectively ν j with θj = θ(d) j = max {θ ij } 1≤i≤d
, to guarantee convergence of the updates.

β-potentials and β-divergences

Assumptions (A) hold for the β-potentials (π β -βπ +β -1)/(β(β -1)) with 0 < β < 1, associated to the so-called β-divergences. Hence, the ROT problem can be solved with the ASA scheme, and the Newton-Raphson steps to update the alternate projections in the POCS technique can be written as follows:

θ ← θ - ((β -1)θ + 1) 1 β-1 1 -p ((β -1)θ + 1) 1 β-1 -1 1 1 , ( 132
) θ ← θ -1 1 ((β -1)θ + 1) 1 β-1 -q 1 ((β -1)θ + 1) 1 β-1 -1 . ( 133
)
Each step can be optimized by computing first the temporary matrix (β -1)θ + 1, then applying an element-wise matrix power of 1/(β -1) -1 to this temporary matrix to obtain a matrix for the denominator, and lastly performing an element-wise matrix multiplication of these two matrices to obtain a matrix for the numerator and thus save one element-wise matrix power. Since ψ is convex and strictly increasing with ψ positive, the convergence of the updates is guaranteed.

Interestingly, the regularizer tends to minus the Burg and Boltzmann-Shannon entropies in the limit β = 0 and β = 1, respectively. Therefore, the β-divergences interpolate between the Itakura-Saito and Kullback-Leibler divergences. We finally remark that the regularizer can also be defined for other values of the parameter β using the same formula, but do not verify assumptions (A) for these values.

p quasi-norms

Considering regularizers -π p with 0 < p < 1, all assumptions (A) are verified except from (A5) since R d×d -⊂ dom ψ = R d×d --. Hence, our primal formulation does not hold here because 0 / ∈ dom ∇ψ. However, it is straightforward to check that our dual formulation for ROT problems with the ASA scheme can still be applied as long as the cost matrix γ does not have null entries so that -γ/λ ∈ dom ∇ψ. Eventually, the Newton-Raphson steps to update the alternate projections in the POCS technique can be written as follows:

θ ← θ + (-θ ) 1 p-1 1 -p 1 p-1 p 1 p-1 (-θ ) 1 p-1 -1 1 1 , ( 134
) θ ← θ + 1 1 (-θ ) 1 p-1 -p 1 p-1 q 1 p-1 1 (-θ ) 1 p-1 -1 . ( 135
)
Each step can be optimized by computing first the temporary matrix -θ , then applying an element-wise matrix power of 1/(p -1) -1 to obtain a matrix for the denominator, and lastly performing an element-wise matrix multiplication of these two matrices to obtain a matrix for the numerator and thus save one element-wise matrix power. Since ψ is convex and strictly increasing with ψ positive everywhere, the convergence of the updates is guaranteed.

p norms

Assumptions (B) hold for the p norms |π| p with 1 < p < +∞, so the ROT problem can be solved with the NASA scheme. For p = 2, the Newton-Raphson steps to update the alternate projections in Dykstra's algorithm can be written as follows:

τ ← τ + sgn(θ -τ 1 ) |θ -τ 1 | 1 p-1 1 -p 1 p-1 p 1 p-1 |θ -τ 1 | 1 p-1 -1 1 , ( 136
) σ ← σ + 1 sgn(θ -τ 1 ) |θ -τ 1 | 1 p-1 -p 1 p-1 q 1 p-1 1 |θ -τ 1 | 1 p-1 -1 . ( 137
)
Denoting θ = θτ 1 or θ = θ -1σ in the respective updates, each step can be optimized by computing first the temporary matrix |θ|, then applying an element-wise matrix power of 1/(p -1) -1 to obtain a matrix for the denominator, and lastly performing an elementwise matrix multiplication of these two matrices and of sgn θ to obtain a matrix for the numerator and thus save one element-wise matrix power as well as several vector replications and matrix subtractions. However, even if ψ is strictly increasing with ψ > 0 on R * , ψ is neither convex nor concave and does not verify the necessary and sufficient condition [START_REF] Rolet | Fast dictionary learning with a smoothed Wasserstein loss[END_REF] for global convergence of the Newton-Raphson method. Moreover, ψ vanishes at 0 for p < 2, and ψ is not differentiable at 0 for p > 2.

Nevertheless, ψ is concave on R -and convex on R + for p < 2, as well as convex on R -and concave on R + for p > 2. It thus divides for a given 1

≤ i ≤ d, respectively 1 ≤ j ≤ d, the real line into at most d + 1 intervals -∞ < θ(1) i ≤ θ(2) i ≤ • • • ≤ θ(d-1) i ≤ θ(d) i < +∞, respectively -∞ < θ(1) j ≤ θ(2) j ≤ • • • ≤ θ(d-1) j ≤ θ(d) j
< +∞, with the values ( θ(k) i ) 1≤k≤d from row i of θ, respectively ( θ(k) j ) 1≤k≤d from column j of θ, sorted in increasing order. The necessary and sufficient condition [START_REF] Rolet | Fast dictionary learning with a smoothed Wasserstein loss[END_REF] is verified on the interior of each of these intervals since we can decompose f (µ i ), respectively g(ν j ), as the sum of an increasing convex and an increasing concave function. Hence, we have global convergence on the interior of the interval that contains the solution. In both cases, we must remove the finite endpoints to ensure differentiability of ψ and positivity of ψ . It is also further possible to prune the last interval from the search. Indeed, we have

d j=1 ψ (θ ij - θ(d) i ) ≤ d j=1 ψ (0) = 0, so that µ i < θi = θ(d) i = max {θ ij } 1≤j≤d . Similarly, we have d i=1 ψ (θ ij - θ(d) j ) ≤ d i=1 ψ (0) = 0, so that ν j < θj = θ(d) j = max {θ ij } 1≤i≤d
. Lastly, we can restrict the first interval with a finite lower bound instead. Indeed, we have

d j=1 ψ (θ ij - θ(1) i + φ (p i /d)) ≥ d j=1 ψ (φ (p i /d)) = p i , so that µ i ≥ θ(1) i -φ (p i /d).
Similarly, we have

d i=1 ψ (θ ij - θ(1) j + φ (q j /d)) ≥ d i=1
ψ (φ (q j /d)) = q j , so that ν j ≥ θ(1) j -φ (q j /d). As a result, we can perform at most d binary searches in parallel to determine within which of the remaining bounded intervals the solutions µ i , respectively ν j , lie. Initialization is then done with the midpoint to guarantee convergence of the updates. A given search thus requires a worst-case logarithmic number of tests, each of which requires a linear number of operations, for a total complexity in O(d 2 log d) instead of O(d 2 ) if no such binary search were needed. Now for p = 2, the regularizer specializes to the Euclidean norm, leading to the squared Euclidean distance as the associated divergence. In addition, the formula for ψ still holds with the convention 0 0 = 1, and ψ is actually constant equal to 1/2. Eventually, the projections can be written in closed form, and we can resort to the analytical algorithm derived in the next example specifically for the Euclidean distance, after doubling the penalty λ to account for the regularizer being halved.

Euclidean Norm and Euclidean Distance

Assumptions (B) hold for half the Euclidean norm π 2 /2 associated to half the squared Euclidean distance. Therefore, the ROT problem can be solved with the NASA scheme, where Dykstra's algorithm can actually be written in closed form. Specifically, the non-negative projection reduces to:

π ← max{0, π} , (138)
and is interleaved with the scaling projections which amount to offsetting the rows and columns of π by an amount such that the rows and columns of π sum up to p and q respectively:

π ← π - 1 d (π 1 -p) 1 , ( 139
) π ← π - 1 d 1 (1 π -q ) . ( 140
)
As a remark, we notice that half the squared Euclidean distance can be seen as a β-divergence using the provided formula for β = 2. However, the β-divergence generated is not of Legendre type because the domain is restricted to R + , whereas it could actually be extended to R so that the regularizer would then be of Legendre type. This is why we fall under assumptions (B) rather than assumptions (A) in this case.

Hellinger Distance

Assumptions (B) hold for the regularizer -(1 -π 2 ) 1 2 akin to a Hellinger distance. Hence, the ROT problem can be solved with the NASA scheme, and the Newton-Raphson steps to update the alternate projections in Dykstra's algorithm can be written as follows:

τ ← τ + (θ -τ 1 ) 1 + (θ -τ 1 ) 2 -1 2 1 -p 1 + (θ -τ 1 ) 2 -3 2 1 , ( 141
) σ ← σ + 1 (θ -1σ ) 1 + (θ -1σ ) 2 -1 2 -q 1 1 + (θ -1σ ) 2 -3 2 . ( 142
)
Denoting θ = θτ 1 or θ = θ -1σ in the respective updates, each step can be optimized by computing first the temporary matrix 1/(1 + θ 2 ), then applying an element-wise matrix square root to this temporary matrix, performing an element-wise matrix multiplication of these two matrices to obtain a matrix for the denominator, and lastly an element-wise matrix multiplication of the temporary matrix with θ to obtain a matrix for the numerator and thus save one element-wise matrix power as well as several vector replications and matrix subtractions. However, even if ψ is strictly increasing with ψ positive everywhere, ψ is neither convex nor concave and does not verify the necessary and sufficient condition [START_REF] Rolet | Fast dictionary learning with a smoothed Wasserstein loss[END_REF] for global convergence of the Newton-Raphson method.

Nevertheless, ψ is convex on R -and concave on R + . It thus divides for a given 1

≤ i ≤ d, respectively 1 ≤ j ≤ d, the real line into at most d + 1 intervals -∞ < θ(1) i ≤ θ(2) i ≤ • • • ≤ θ(d-1) i ≤ θ(d) i < +∞, respectively -∞ < θ(1) j ≤ θ(2) j ≤ • • • ≤ θ(d-1) j ≤ θ(d) j < +∞, with
the values ( θ(k) i ) 1≤k≤d from row i of θ, respectively ( θ(k) j ) 1≤k≤d from column j of θ, sorted in increasing order. On each of these intervals, the necessary and sufficient condition ( 40) is verified since we can decompose f (µ i ), respectively g(ν j ), as the sum of an increasing convex and an increasing concave function. Hence, we have global convergence on the interval that contains the solution. It is further possible to prune the last interval from the search. Indeed, we have

d j=1 ψ (θ ij - θ(d) i ) ≤ d j=1 ψ (0) = 0, so that µ i < θi = θ(d) i = max {θ ij } 1≤j≤d . Similarly, we have d i=1 ψ (θ ij - θ(d) j ) ≤ d i=1 ψ (0) = 0, so that ν j < θj = θ(d) j = max {θ ij } 1≤i≤d
. Lastly, we can restrict the first interval with a finite lower bound instead. Indeed, we have

d j=1 ψ (θ ij - θ(1) i + φ (p i /d)) ≥ d j=1 ψ (φ (p i /d)) = p i , so that µ i ≥ θ(1) i -φ (p i /d). Similarly, we have d i=1 ψ (θ ij - θ(1) j + φ (q j /d)) ≥ d i=1 ψ (φ (q j /d)) = q j ,
so that ν j ≥ θ(1) j -φ (q j /d). As a result, we can perform d binary searches in parallel to determine within which of the remaining intervals the solutions µ i , respectively ν j , lie. Initialization is then done with the midpoint to guarantee convergence of the updates. A given search requires a worst-case logarithmic number of tests, each of which requires a linear number of operations, for a total complexity in O(d 2 log d) instead of O(d 2 ) if no such binary search were needed.

Quadratic Forms and Mahalanobis Distances

Assumptions (B) hold for the quadratic forms (1/2) vec(π) P vec(π) with positive-definite matrix P ∈ R d 2 ×d 2 , associated to the Mahalanobis distances, so the ROT problem can be solved with the NASA scheme. For a diagonal matrix P, the regularizer is separable and the Newton-Raphson steps to update the alternate projections in Dykstra's algorithm are similar to that for the Euclidean distance with appropriate weights. For a non-diagonal matrix P, however, the regularizer is not separable anymore and we must resort to the generic NASA scheme.

In this general case, the scaling projections amount to convex quadratic programs with linear equality constraints. They can be solved using classical techniques such as the range-space and null-space approaches, Krylov subspace methods or active set strategies. The non-negative projection reduces to a convex quadratic program with a linear inequality constraint. It can be solved elegantly with an iterative algorithm for non-negative quadratic programming proposed by [START_REF] Sha | Multiplicative updates for nonnegative quadratic programming[END_REF] using multiplicative updates with a complexity in O(d 4 ). All in all, we recommend using a sparse matrix P with a block-diagonal structure and an order of magnitude of d 2 non-null entries, so as to obtain a quadratic instead of quartic empirical complexity.

Experimental Results

In this section, we present the results of our methods on different experiments. We first design an synthetic test to showcase the behavior of different regularizers and penalties on the output solutions or computational times (Section 6.1). We then consider a pattern recognition application to audio scene classification on a real-world dataset (Section 6.2).

Synthetic Data

We start by visualizing the effects of different regularizers φ and varying penalties λ on synthetic data. For the input distributions, we discretize and normalize continuous densities on a uniform grid (x i ) 1≤i≤d of [0, 1] with dimension d = 256. We use for p a univariate normal with mean 0.5 and variance 0.2, and for q a mixture of two normals with equal weights, respective means 0.25 and 0.75, and same variance 0.1. We set the cost matrix γ as the squared Euclidean distance

γ ij = (x i -x j )
2 on the grid. The input distributions (bottom left and top right), cost matrix (top left) and unique earth mover's plan (bottom right) computed for classical OT using the solver of [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF] with standard settings, are shown in Figure 3.

We test all separable regularizers φ introduced in Section 5. Because these regularizers have different ranges in the sensible values of the rot mover's plans π , we manually tune the penalties λ so that they feature similar amounts of regularization. For ease of comparison, we set λ = λ λ , with λ constant for each φ, and λ varying similarly for all φ. The limit case when λ tends to infinity is simply obtained by setting γ/λ = 0 in the algorithms, except from p quasi-norms for which we use λ = 10 10 . The null values of γ are also fixed to 10 -12 for p quasi-norms. We do Figure 3: Earth mover's plan π for the cost matrix γ and input distributions p, q.

not limit the number of iterations in the different algorithms, and use a small tolerance of 10 -8 for convergence with the ∞ norm on the marginal difference checked after each iteration as a termination criterion.

The rot mover's plans obtained for ROT for d = 256 with the different regularizers and penalties are visualized in Figure 4. We first observe that all rot mover's plans converge to the earth mover's plan for low values of the penalty as shown theoretically in Property 9. Nevertheless, the rot mover's plans exhibit different shapes depending on the regularizers for intermediary and large values of the penalty. In the limit when the penalty grows to infinity, we obtain the transport plan with minimal Bregman information as shown theoretically in Property 7. In particular, this leads to pq with an ellipsoidal shape for BSKL (Boltzmann-Shannon entropy and Kullback-Leibler divergence), meaning that the mass is relatively spread among neighbor bins. The same pattern is observed for FDLOG (Fermi-Dirac entropy and logistic loss function), which can be explained in this synthetic example by the rot mover's plans having low values and the two regularizers being equivalent up to a constant in the neighborhood of zero. The profile gets more rectangular for BIS (Burg entropy and Itakura-Saito divergence), implying that the mass is even more spread across the different bins. Using an intermediary value β = 0.5 in BETA (β-potentials and β-divergences) allows the interpolation between these two limits of a rectangle for β = 1 and an ellipsoid for β = 0, so that the parameter β actually helps to control the spread of mass in the regularization. We observe similar results for LPQN ( p quasi-norms) with an ellipsoid for p = 0.9, a rectangle for p = 0.1, and a shape in between for p = 0.5. When the power parameter further increases in LPN ( p norms), we obtain new shapes that feature less spread of mass. These shapes for p = 1.1 and p = 1.5 now interpolate up to a lozenge for p = 2 in EUC (Euclidean norm and Euclidean distance), so that the parameter p also provides control on the spread of mass. A similar diamond profile is obtained for HELL (Hellinger distance), which is due again to the rot mover's plans having low values and the two regularizers being equivalent up to a constant in the neighborhood of zero. Lastly, we remark that varying the penalty between the two extremes allows a smooth interpolation of the earth mover's plan and optimal plan with minimal Bregman information, while keeping similar shapes and effects in terms of spreading of mass.

We next report in Table 4 the computational times required to reach convergence for the different regularizers and penalties. As a stopping criterion, we use the relative variation with tolerance 10 -2 in 2 norm for the main loop of alternate Bregman projections, and the absolute variation with tolerance 10 -5 in 2 norm for the auxiliary loops of the Newton-Raphson method. We use the same synthetic data as above but also vary the dimension d to assess its influence on speed. As already observed specifically for Sinkhorn distances [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF], computing ROT distances is faster for important regularization with larger values of λ. The regularizers under assumptions (A) do not require the extra projections onto the non-negative orthant, and thus intuitively require less computational effort than the ones that verify assumptions (B). In addition, we notice that when the projections have closed-form expressions, the algorithms are also faster.

The results further illustrate the influence of the data dimension d and the difference between ROT and classical OT performances. For a low dimension d, the RMD is competitive with EMD in his historical implementation EMD OLD [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF]. The super-cubic complexity of the EMD with EMD OLD becomes prohibitive as the data dimension increases in contrast to the RMD which scales better. It should nevertheless be underlined that for reasonable dimensions, fast computation of the EMD can be obtained with a more recent, optimized implementation of the network simplex solver EMD NEW [START_REF] Bonneel | Displacement interpolation using lagrangian mass transport[END_REF]. For higher dimensions, the super-cubic complexity makes EMD NEW less attractive, though it stays competitive with the RMD under a dimension d = 4096. As a consequence, a numerical alternative to our algorithms for solving ROT problems with reasonable dimensions is to rely on conditional gradient methods similar to [START_REF] Ferradans | Regularized discrete optimal transport[END_REF]. Indeed, such methods imply the iterative resolution of linearized ROT problems, that can be reformulated as EMD problems and therefore be solved with the fast network simplex approach [START_REF] Bonneel | Displacement interpolation using lagrangian mass transport[END_REF]. Lastly, for a fair interpretation of the above timing results, we must mention that the two EMD schemes tested were run under MATLAB from native C/C++ implementations 12 via compiled MEX files 34 . Hence, these EMD codes are quite optimized in comparison to our pure MATLAB prototype codes5 for the RMD. It is thus plausible that optimized C/C++ implementations of our algorithms would be even more competitive in this context.

Audio Classification

We now assess our methods in the context of audio classification, and specifically address the task of acoustic scene classification where the goal is to assign a test recording to one of predefined classes that characterizes the environment in which it was captured. We consider the framework of the DCASE 2016 IEEE AASP challenge with the TUT Acoustic Scenes 2016 database [START_REF] Mesaros | TUT database for acoustic scene classification and sound event detection[END_REF]. A baseline system is also provided with the database for comparison. This system is based on Mel-frequency cepstral coefficient (MFCC) timbral features with Gaussian mixture model (GMM) classification. One GMM with diagonal covariance matrix is learned per class by expectation-maximization (EM), after concatenating and normalizing in mean and variance the extracted MFCCs from the training segments in that class. A test file is assigned to the class whose trained GMM leads to maximum likelihood for the extracted MFCCs for that file, where the MFCCs are considered as independent samples and normalized with the learned mean and variance for the respective classes. The baseline system is ran with its default parameters: 40 ms frame size, 20 ms hop size, 60-dimensional MFCCs comprising 20 static (including energy) plus 20 delta and 20 acceleration coefficients extracted with standard settings in RASTAMAT, 16 GMM components learned with standard settings in VOICEBOX.

Since MFCCs potentially take negative values, OT tools cannot be applied directly to this kind of features. Therefore, the common approach is to compute OT appropriately on GMMs estimated from MFCCs instead. Our proposed system follows this principle, and is implemented in the very same pipeline as the baseline for a fair comparison, with the following differences. One GMM is learned by EM for each training segment instead of class. Any normalization on the MFCCs per class is thus removed. Since less components are typically required to model one segment compared to one class, the spurious GMM components are further discarded as postprocessing by keeping only those with weight and variances all greater than 10 -2 . Instead of applying a GMM classifier, all individual models are exploited to train a support vector machine (SVM) classifier. An exponential kernel for the SVM is designed by introducing a distance between two mixtures P, Q based on the RMD as follows: κ(P, Q) = exp(-d γ,λ,φ (ω, υ)/τ ) , (

where the exponential decay rate τ > 0 is a kernel parameter, and ω, υ ∈ Σ d are the respective weights of the d = 16 (or less) components for the two GMMs P, Q. The cost matrix γ ∈ R d×d + depends on P, Q and is the square root of a symmetrized Kullback-Leibler divergence, called the Jeffrey divergence, between the pairwise Gaussian components:

γ ij = 1 4 l k=1 (σ 2 ik -ς 2 jk ) 2 + (σ 2 ik + ς 2 jk )(µ ik -ν jk ) 2 σ 2 ik ς 2 jk , ( 144
)
where µ i and σ 2 i , respectively ν j and ς 2 j , are the means and variances of the l = 60 MFCC features for component i in the first mixture P , respectively component j in the second mixture Q. The SVM classifier is implemented with standard settings in LIBSVM, and requires an additional soft-margin parameter C > 0 to be tuned. Notice that, even if the kernel is not positivedefinite, LIBSVM is still able to provide a relevant classification by guaranteeing convergence to a stationary point [START_REF] Alabdulmohsin | Support vector machines with indefinite kernels[END_REF][START_REF] Haasdonk | Feature space interpretation of SVMs with indefinite kernels[END_REF][START_REF] Lin | A study on sigmoid kernels for SVM and the training of non-PSD kernels by SMO-type methods[END_REF]. All separable regularizers φ from Section 5 with different penalties λ > 0 are tested for the RMD in comparison to the EMD. The two distances between p, q and q, p with cost matrix γ transposed are computed and averaged, so as to remove any asymmetry due to practical issues. The number of iterations is limited to 100 for the main loop of the algorithm and to 10 for the auxiliary loops of the Newton-Raphson method, and the tolerance is set to 10 -6 in all loops for convergence with the ∞ norm on the marginal difference checked after each iteration as a termination criterion. The parameters τ, C ∈ 10 {-1,+0,+1,+2} and penalty λ ∈ Λ, where Λ is a manually chosen set of four successive powers of ten depending on the range of the regularizer φ, are tuned automatically by cross-validation.

The obtained results on this experiment in terms of accuracy per system are reported in Table 5. The optimal penalties λ ∈ Λ selected by cross-validation for each regularizer φ are also included, while the optimal parameters τ, C are not displayed since they actually all equal 10 +1 independently of the kernel used. We first notice that the proposed system SVM (support vector machine classifier) consistently outperforms the baseline system GMM (Gaussian mixture model classifier). This proves the benefits of incorporating individual information per sound via an SVM rather than exploiting global information per class with a GMM. This further demonstrates the relevance of OT and more general ROT problems for the design of kernels between GMMs in the SVM pipeline. We also notice that RMD (rot mover's distance kernel) is at least competitive with EMD (earth mover's distance kernel) for all proposed regularizers, except from EUC which does not perform as well. This might be a consequence of the regularization profile for EUC, or equivalently LPN with p = 2, which does not spread enough mass across similar bins, implying a lack of robustness to slight variations in the means and variances of the GMM components. Reducing the power parameter in LPN brings back to a competitive system with EMD for p = 1.1, and even a better trade-off with improved accuracy for p = 1.5. We obtain similar results for LPQN with p = 0.9 and p = 0.5, with now the best compromise for the lowest power value p = 0.1 which clearly outperforms EMD. As a remark, the accuracy for LPN and LPQN is not unimodal with respect to p which controls the spread of mass in the regularization. We suspect this is because the performance is a function of both the spread of mass and the amount of regularization, whose coupling allows for similar compromises in terms of results within different regimes of use. Concerning BETA now, we observe that the existing Sinkhorn-Knopp algorithm BSKL for β = 1 does not improve the accuracy compared to EMD. 1.00 10 {-2,-1,+0,+1} 10 +0 81.3% -BETA 0.50 10 {-3,-2,-1,+0} 10 -2 81.5% BIS BETA 0.00 10 {-4,-3,-2,-1} 10 -2 81.3% -LPQN 0.10 10 {-2,-1,+0,+1} 10 -1 82.1% -LPQN 0.50 10 {-2,-1,+0,+1} 10 -1 81.3% -LPQN 0.90 10 {-2,-1,+0,+1} Increasing the spread of mass with β = 0 in BIS is even worse. The best performance is obtained with a range in between for β = 0.5, which slightly improves results over EMD. Using LOG here slightly degrades the performance compared to EMD and BSKL. Interestingly, the overall best accuracy on this application is obtained for HELL which beats all other systems, including EUC, by a safe margin. In contrast to the experiment on synthetic data with dimension 256 presented in Section 6.1, where both BSKL and LOG, respectively EUC and HELL, behave similarly due to equivalence up to a constant for low values in the transport plans, the range of the transport plans here is much higher since the dimension of the input distributions is at most 16 (typically less than 10). This raises the importance of choosing a good regularizer depending on the actual task and its inherent design criteria such as the data dimension.

Conclusion

In this paper, we formulated a unified framework for smooth convex regularization of discrete OT problems. We also derived some algorithmic methods to solve such ROT problems, and detailed their specificities for classical regularizers and associated divergences from the literature. We finally designed a synthetic experiment to illustrate our proposed methods, and proved the relevance of ROT problems and the RMD on a real-world application to audio scene classification. The obtained results are encouraging for further development of the present work, and we now discuss some interesting perspectives for future investigation. Firstly, we want to assess the effect of other regularizers on the solutions, notably when adding an affine term. From a geometrical viewpoint, such a transformation is equivalent to simply translating the cost matrix, with no effect on the Bregman divergence itself. For a given regularizer, we could therefore parametrize a whole family of interpolating regularizers, and tune the translation parameter according to the application. In particular, a recent work developed independently of ours makes use of Tsallis entropies to regularize OT problems with ad hoc solvers [START_REF] Muzellec | Tsallis regularized optimal transport and ecological inference[END_REF]. These regularizers could be integrated readily to our more general framework based on alternate Bregman projections, since Tsallis entropies are equivalent to β-potentials and p (quasi)-norms up to an affine term.

In another direction, we would like to extend some theoretical results that hold for the Boltzmann-Shannon entropy and associated Kullback-Leibler divergence. Specifically, it is known that the related rot mover's plan converges in norm to the earth mover's plan with an exponential rate as the penalty decreases [START_REF] Cominetti | Asymptotic analysis of the exponential penalty trajectory in linear programming[END_REF]. It is not straightforward, however, to generalize this to other regularizers and divergences. In addition, it would be worth elucidating some technical restrictions under which metric properties such as the triangular inequality can be proved similarly to Sinkhorn distances.

We also plan to study other pattern recognition tasks in text, image and audio signal processing. Intuitive possibilities include retrieval and classification for various kinds of data modeled via histograms of features or GMMs. Among potential approaches, this can be addressed by exploiting the RMD either directly in a nearest-neighbor search, or in the design of kernels for an SVM as done here for acoustic scenes. For such tasks, it would be relevant to provide insight into the choice of a good regularizer for the actual problem, or develop methods for automatic tuning of regularization parameters, and for learning the cost matrix from the data as can be done for the EMD [START_REF] Cuturi | Ground metric learning[END_REF]. Even if we mostly focused on separable regularizers, it would be relevant to further use the quadratic forms associated to Mahalanobis distances in certain applications, and maybe propose a parametric learning scheme for the quadratic regularizer from the data.

Lastly, a more prospective idea is to use the RMD instead of Sinkhorn distances in the recent works built on the entropic regularization mentioned in Section 1. We also think that variational ROT problems could be formulated for statistical inference, notably parameter estimation in finite mixture models by minimizing loss functions based on the RMD [START_REF] Dessein | Parameter estimation in finite mixture models by regularized optimal transport: A unified framework for hard and soft clustering[END_REF]. This would leverage new applications of our ROT framework for more general machine learning problems. Such developments are yet involved and require some theoretical effort before reaching enough maturity to address practical setups.

Figure 4 :

 4 Figure 4: Rot mover's plans π for different regularizers φ and penalties λ = λ λ .

Table 4 :

 4 10 -2 10 -1 10 +0 10 -2 10 -1 10 +0 10 -2 10 -1 10 +0 Computational times in seconds required to reach convergence for different regularizers φ and penalties λ = λ λ , with varying dimensions d.

	10 +0

  The data set consists of audio recordings at 44.1 kHz sampling rate and 24-bit resolution. The metadata contains ground-truth annotations on the type of acoustic scene for all files, with a total of 15 classes: home, office, library, café/restaurant, grocery store, city center, residential area, park, forest path, beach, car, train, bus, tram, metro station. The audio material is cut into 30second segments, and is split into two subsets of 75%-25% containing respectively 78-26 segments per class for development and evaluation, resulting in a total of 1170-390 files for training and testing. A 4-fold cross-validation setup is given with the training set. The classification accuracy, that is, the number of correctly classified segments among the total number of segments, is used as a score to evaluate systems.

Table 5 :

 5 Results of the experiment on audio classification.

	10 +0	81.0%

http://robotics.stanford.edu/~rubner/emd/default.htm

http://liris.cnrs.fr/~nbonneel/FastTransport/

https://github.com/francopestilli/life/tree/master/external/emd

https://arolet.github.io/code/

https://www.math.u-bordeaux.fr/~npapadak/GOTMI/codes.php

Acknowledgments This study has been carried out with financial support from the French State, managed by the French National Research Agency (ANR) in the frame of the "Investments for the future" Program IdEx Bordeaux (ANR-10-IDEX-03-02), Cluster of excellence CPU and the GOTMI project (ANR-16-CE33-0010-01). The authors would like to thank Annamaria Mesaros for her kind help with the evaluation on the DCASE 2016 IEEE AASP challenge, Charles-Alban Deledalle for his valuable advice on the use of the computing platform PlaFRIM, and Marco Cuturi for the insightful discussions about this work.