Mokhtar Z Alaya
email: mokhtarzahdi.alaya@gmail.com

Litis Ea4108

Maxime Bérar
email: maxime.berar@univ-rouen.fr

Gilles Gasso
email: gilles.gasso@insa-rouen.fr

Alain Rakotomamonjy

Screening Sinkhorn Algorithm for Regularized Optimal Transport

We introduce in this paper a novel strategy for efficiently approximating the Sinkhorn distance between two discrete measures. After identifying neglectable components of the dual solution of the regularized Sinkhorn problem, we propose to screen those components by directly setting them at that value before entering the Sinkhorn problem. This allows us to solve a smaller Sinkhorn problem while ensuring approximation with provable guarantees. More formally, the approach is based on a new formulation of dual of Sinkhorn divergence problem and on the KKT optimality conditions of this problem, which enable identification of dual components to be screened. This new analysis leads to the SCREENKHORN algorithm. We illustrate the efficiency of SCREENKHORN on complex tasks such as dimensionality reduction and domain adaptation involving regularized optimal transport.

Introduction

Computing optimal transport (OT) distances between pairs of probability measures or histograms, such as the earth mover's distance [START_REF] Werman | A distance metric for multidimensional histograms[END_REF][START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF] and Monge-Kantorovich or Wasserstein distance [START_REF] Villani | Optimal Transport: Old and New[END_REF], are currently generating an increasing attraction in different machine learning tasks [START_REF] Solomon | Wasserstein propagation for semi-supervised learning[END_REF][START_REF] Kusner | From word embeddings to document distances[END_REF][START_REF] Arjovsky | Wasserstein generative adversarial networks[END_REF][START_REF] Ho | Multilevel clustering via Wasserstein means[END_REF], statistics [START_REF] Frogner | Learning with a Wasserstein loss[END_REF][START_REF] Panaretos | Amplitude and phase variation of point processes[END_REF][START_REF] Ebert | Construction of non-asymptotic confidence sets in 2-Wasserstein space[END_REF][START_REF] Bigot | Geodesic PCA in the Wasserstein space by convex PCA[END_REF][START_REF] Flamary | Wasserstein discriminant analysis[END_REF], and computer vision [START_REF] Bonneel | Displacement interpolation using Lagrangian mass transport[END_REF][START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF][START_REF] Solomon | Convolutional Wasserstein distances: Efficient optimal transportation on geometric domains[END_REF], among other applications [START_REF] Kolouri | Optimal mass transport: Signal processing and machine-learning applications[END_REF][START_REF] Peyré | Computational optimal transport[END_REF]. In many of these problems, OT exploits the geometric features of the objects at hand in the underlying spaces to be leveraged in comparing probability measures. This effectively leads to improved performance of methods that are oblivious to the geometry, for example the chi-squared distances or the Kullback-Leibler divergence. Unfortunately, this advantage comes at the price of an enormous computational cost of solving the OT problem, that can be prohibitive in large scale applications. For instance, the OT between two histograms with supports of equal size n can be formulated as a linear programming problem that requires generally super O(n 2.5 ) [START_REF] Lee | Path finding methods for linear programming: Solving linear programs in Õ(vrank) iterations and faster algorithms for maximum flow[END_REF] arithmetic operations, which is problematic when n becomes larger.

A remedy to the heavy computation burden of OT lies in a prevalent approach referred to as regularized OT [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] and operates by adding an entropic regularization penalty to the original problem.

33rd Conference on Neural Information Processing Systems (NeurIPS 2019), Vancouver, Canada.

Such a regularization guarantees a unique solution, since the objective function is strongly convex, and a greater computational stability. More importantly, this regularized OT can be solved efficiently with celebrated matrix scaling algorithms, such as Sinkhorn's fixed point iteration method [START_REF] Sinkhorn | Diagonal equivalence to matrices with prescribed row and column sums[END_REF][START_REF] Knight | The Sinkhorn-Knopp algorithm: Convergence and applications[END_REF][START_REF] Kalantari | On the complexity of general matrix scaling and entropy minimization via the ras algorithm[END_REF].

Several works have considered further improvements in the resolution of this regularized OT problem. A greedy version of Sinkhorn algorithm, called Greenkhorn Altschuler et al. (2017), allows to select and update columns and rows that most violate the polytope constraints. Another approach based on low-rank approximation of the cost matrix using the Nyström method induces the Nys-Sink algorithm [START_REF] Altschuler | Massively scalable Sinkhorn distances via the Nyström method[END_REF]. Other classical optimization algorithms have been considered for approximating the OT, for instance accelerated gradient descent [START_REF] Xie | A fast proximal point method for computing Wasserstein distance[END_REF][START_REF] Dvurechensky | Computational optimal transport: Complexity by accelerated gradient descent is better than by Sinkhorn's algorithm[END_REF][START_REF] Lin | On efficient optimal transport: An analysis of greedy and accelerated mirror descent algorithms[END_REF], quasi-Newton methods [START_REF] Blondel | Smooth and sparse optimal transport[END_REF][START_REF] Cuturi | A smoothed dual approach for variational Wasserstein problems[END_REF] and stochastic gradient descent [START_REF] Genevay | Stochastic optimization for large-scale optimal transport[END_REF][START_REF] Abid | Stochastic algorithms for entropy-regularized optimal transport problems[END_REF].

In this paper, we propose a novel technique for accelerating the Sinkhorn algorithm when computing regularized OT distance between discrete measures. Our idea is strongly related to a screening strategy when solving a Lasso problem in sparse supervised learning [START_REF] Ghaoui | Safe feature elimination in sparse supervised learning[END_REF]. Based on the fact that a transport plan resulting from an OT problem is sparse or presents a large number of neglectable values [START_REF] Blondel | Smooth and sparse optimal transport[END_REF], our objective is to identify the dual variables of an approximate Sinkhorn problem, that are smaller than a predefined threshold, and thus that can be safely removed before optimization while not altering too much the solution of the problem. Within this global context, our contributions are the following:

• From a methodological point of view, we propose a new formulation of the dual of the Sinkhorn divergence problem by imposing variables to be larger than a threshold. This formulation allows us to introduce sufficient conditions, computable beforehand, for a variable to strictly satisfy its constraint, leading then to a "screened" version of the dual of Sinkhorn divergence. • We provide some theoretical analysis of the solution of the "screened" Sinkhorn divergence, showing that its objective value and the marginal constraint satisfaction are properly controlled as the number of screened variables decreases. • From an algorithmic standpoint, we use a constrained L-BFGS-B algorithm [START_REF] Nocedal | Updating quasi-newton matrices with limited storage[END_REF][START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF] but provide a careful analysis of the lower and upper bounds of the dual variables, resulting in a well-posed and efficient algorithm denoted as SCREENKHORN.

• Our empirical analysis depicts how the approach behaves in a simple Sinkhorn divergence computation context. When considered in complex machine learning pipelines, we show that SCREENKHORN can lead to strong gain in efficiency while not compromising on accuracy.

The remainder of the paper is organized as follow. In Section 2 we briefly review the basic setup of regularized discrete OT. Section 3 contains our main contribution, that is, the SCREENKHORN algorithm. Section 4 is devoted to theoretical guarantees for marginal violations of SCREENKHORN.

In Section 5 we present numerical results for the proposed algorithm, compared with the state-of-art Sinkhorn algorithm as implemented in [START_REF] Flamary | POT: Python optimal transport library[END_REF]. The proofs of theoretical results are postponed to the supplementary material as well as additional empirical results.

Notation. For any positive matrix T ∈ R n×m , we define its entropy as

H(T ) = -i,j T ij log(T ij ). Let r(T ) = T 1 m ∈ R n and c(T ) = T 1 n ∈ R m denote
the rows and columns sums of T respectively. The coordinates r i (T ) and c j (T ) denote the i-th row sum and the j-th column sum of T , respectively. The scalar product between two matrices denotes the usual inner product, that is

T, W = tr(T W ) = i,j T ij W ij ,
where T is the transpose of T . We write 1 (resp. 0) the vector having all coordinates equal to one (resp. zero). ∆(w) denotes the diag operator, such that if

w ∈ R n , then ∆(w) = diag(w 1 , . . . , w n ) ∈ R n×n . For a set of indices L = {i 1 , . . . , i k } ⊆ {1, . . . , n} satisfying i 1 < • • • < i k ,

Regularized discrete OT

We briefly expose in this section the setup of OT between two discrete measures. We then consider the case when those distributions are only available through a finite number of samples, that is

µ = n i=1 µ i δ xi ∈ Σ n and ν = m j=1 ν i δ yj ∈ Σ m
, where Σ n is the probability simplex with n bins, namely the set of probability vectors in R n + , i.e., Σ n = {w ∈ R n + :

n i=1 w i = 1}. We denote their probabilistic couplings set as Π(µ, ν) = {P ∈ R n×m + , P 1 m = µ, P 1 n = ν}.
Sinkhorn divergence. Computing OT distance between the two discrete measures µ and ν amounts to solving a linear problem [START_REF] Kantorovich | On the transfer of masses (in russian)[END_REF] given by S(µ, ν) = min P ∈Π(µ,ν) C, P , where P = (P ij ) ∈ R n×m is called the transportation plan, namely each entry P ij represents the fraction of mass moving from x i to y j , and C = (C ij ) ∈ R n×m is a cost matrix comprised of nonnegative elements and related to the energy needed to move a probability mass from x i to y j . The entropic regularization of OT distances [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] relies on the addition of a penalty term as follows:

S η (µ, ν) = min

P ∈Π(µ,ν) { C, P -ηH(P )}, (1)
where η > 0 is a regularization parameter. We refer to S η (µ, ν) as the Sinkhorn divergence [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF].

Dual of Sinkhorn divergence. Below we provide the derivation of the dual problem for the regularized OT problem (1). Towards this end, we begin with writing its Lagrangian dual function:

L (P, w, z) = C, P + η log P, P + w, P 1 m -µ + z, P 1 n -ν .
The dual of Sinkhorn divergence can be derived by solving min P ∈R n×m + L (P, w, z). It is easy to check that objective function P → L (P, w, z) is strongly convex and differentiable. Hence, one can solve the latter minimum by setting ∇ P L (P, w, z) to 0 n×m . Therefore, we get P ij = exp -1 η (w i + z j + C ij ) -1 , for all i = 1, . . . , n and j = 1, . . . , m. Plugging this solution, and setting the change of variables u = -w/η -1/2 and v = -z/η -1/2, the dual problem is given by

min u∈R n ,v∈R m Ψ(u, v) := 1 n B(u, v)1 m -u, µ -v, ν , (2)
where B(u, v) := ∆(e u )K∆(e v ) and K := e -C/η stands for the Gibbs kernel associated to the cost matrix C. We refer to problem (2) as the dual of Sinkhorn divergence. Then, the optimal solution P of the primal problem (1) takes the form P = ∆(e u )K∆(e v ) where the couple (u , v ) satisfies:

(u , v ) = argmin u∈R n ,v∈R m {Ψ(u, v)}.
Note that the matrices ∆(e u ) and ∆(e v ) are unique up to a constant factor [START_REF] Sinkhorn | Diagonal equivalence to matrices with prescribed row and column sums[END_REF]. Moreover, P can be solved efficiently by iterative Bregman projections [START_REF] Benamou | Iterative bregman projections for regularized transportation problems[END_REF] referred to as Sinkhorn iterations, and the method is referred to as SINKHORN algorithm which, recently, has been proven to achieve a near-O(n 2 ) complexity [START_REF] Altschuler | Near-linear time approximation algorithms for optimal transport via Sinkhorn iteration[END_REF]. Motivation. The key idea of our approach is motivated by the so-called static screening test [START_REF] Ghaoui | Safe feature elimination in sparse supervised learning[END_REF] in supervised learning, which is a method able to safely identify inactive features, i.e., features that have zero components in the solution vector. Then, these inactive features can be removed from the optimization problem to reduce its scale. Before diving into detailed algorithmic analysis, let us present a brief illustration of how we adapt static screening test to the dual of Sinkhorn divergence. Towards this end, we define the convex set C r α ⊆ R r , for r ∈ N and α > 0, by C r α = {w ∈ R r : e wi ≥ α}. In Figure 1, we plot (e u , e v ) where (u , v ) is the pair solution of the dual of Sinkhorn divergence (2) in the particular case of: n = m = 500, η = 1, µ = ν = 1 n 1 n , x i ∼ N ((0, 0) , ( 1 0 0 1 )), y j ∼ N ((3, 3) , 1 -0.8 -0.8 1

) and the cost matrix C corresponds to the pairwise euclidean distance, i.e., C ij = x i -y j 2 . We also plot two lines corresponding to e u ≡ α u and e v ≡ α v for some α u > 0 and α v > 0, choosing randomly and playing the role of thresholds to select indices to be discarded. If we are able to identify these indices before solving the problem, they can be fixed at the thresholds and removed then from the optimization procedure yielding an approximate solution.

Static screening test. Based on this idea, we define a so-called approximate dual of Sinkhorn divergence

min u∈C n ε κ ,v∈C m εκ Ψ κ (u, v) := 1 n B(u, v)1 m -κu, µ - v κ , ν , (3)
which is simply a dual of Sinkhorn divergence with lower-bounded variables, where the bounds are α u = εκ -1 and α v = εκ with ε > 0 and κ > 0 being fixed numeric constants which values will be clear later. The new formulation (3) has the form of (κµ, ν/κ)-scaling problem under constraints on the variables u and v. Those constraints make the problem significantly different from the standard scaling-problems [START_REF] Kalantari | On the complexity of nonnegative-matrix scaling[END_REF]. We further emphasize that κ plays a key role in our screening strategy. Indeed, without κ, e u and e v can have inversely related scale that may lead in, for instance e u being too large and e v being too small, situation in which the screening test would apply only to coefficients of e u or e v and not for both of them. Moreover, it is clear that the approximate dual of Sinkhorn divergence coincides with the dual of Sinkhorn divergence (2) when ε = 0 and κ = 1. Intuitively, our hope is to gain efficiency in solving problem (3) compared to the original one in Equation ( 2) by avoiding optimization of variables smaller than the threshold and by identifying those that make the constraints active. More formally, the core of the static screening test aims at locating two subsets of indices (I, J) in {1, . . . , n} × {1, . . . , m} satisfying: e ui > α u , and e vj > α v , for all (i, j) ∈ I × J and e u i = α u , and e v j = α v , for all (i , j

) ∈ I × J , namely (u, v) ∈ C n αu × C m αv .
The following key result states sufficient conditions for identifying variables in I and J . Lemma 1. Let (u * , v * ) be an optimal solution of problem (3). Define

I ε,κ = i = 1, . . . , n : µ i ≥ ε 2 κ r i (K) , J ε,κ = j = 1, . . . , m : ν j ≥ κε 2 c j (K) (4)
Then one has e u * i = εκ -1 and e v * j = εκ for all i ∈ I ε,κ and j ∈ J ε,κ .

Proof of Lemma 1 is postponed to the supplementary material. It is worth to note that first order optimality conditions applied to (u * , v * ) ensure that if e u * i > εκ -1 then e u * i (Ke v * ) i = κµ i and if e v * j > εκ then e v * j (K e u * ) j = κ -1 ν j , that correspond to the Sinkhorn marginal conditions [START_REF] Peyré | Computational optimal transport[END_REF] up to the scaling factor κ.

Screening with a fixed number budget of points. The approximate dual of Sinkhorn divergence is defined with respect to ε and κ. As those parameters are difficult to interpret, we exhibit their relations with a fixed number budget of points from the supports of µ and ν. In the sequel, we denote by n b ∈ {1, . . . , n} and m b ∈ {1, . . . , m} the number of points that are going to be optimized in problem (3), i.e, the points we cannot guarantee that e u * i = εκ -1 and e v * j = εκ .

Let us define ξ ∈ R n and ζ ∈ R m to be the ordered decreasing vectors of µ r(K) and ν c(K) respectively, that is

ξ 1 ≥ ξ 2 ≥ • • • ≥ ξ n and ζ 1 ≥ ζ 2 ≥ • • • ≥ ζ m .
To keep only n b -budget and m b -budget of points, the parameters κ and ε satisfy We are now in position to formulate the optimization problem related to the screened dual of Sinkhorn. Indeed, using the above analyses, any solution (u * , v * ) of problem (3) satisfies e u * i ≥ εκ -1 and e v * j ≥ εκ for all (i, j) ∈ (I ε,κ × J ε,κ ), and e u * i = εκ -1 and e v * j = εκ for all (i, j) ∈ (I ε,κ × J ε,κ ). Hence, we can restrict the problem (3) to variables in I ε,κ and J ε,κ . This boils down to restricting the constraints feasibility C n ε κ ∩ C m εκ to the screened domain defined by U sc ∩ V sc ,

ε 2 κ -1 = ξ n b and ε 2 κ = ζ m b . Hence ε = (ξ n b ζ m b ) 1/4 and κ = ζ m b ξ n b . (5)
U sc = {u ∈ R n b : e u Iε,κ ε κ 1 n b } and V sc = {v ∈ R m b : e v Jε,κ εκ1 m b }
where the vector comparison has to be understood elementwise. And, by replacing in Equation ( 3), the variables belonging to (I ε,κ × J ε,κ ) by εκ -1 and εκ, we derive the screened dual of Sinkhorn divergence problem as

min u∈Usc,v∈Vsc {Ψ ε,κ (u, v)} (6)
where

Ψ ε,κ (u, v) = (e u Iε,κ ) K (Iε,κ,Jε,κ) e v Jε,κ + εκ(e u Iε,κ ) K (Iε,κ,J ε,κ ) 1 m b + εκ -1 1 n b K (I ε,κ ,Jε,κ) e v Jε,κ -κµ Iε,κ u Iε,κ -κ -1 ν Jε,κ v Jε,κ + Ξ with Ξ = ε 2 i∈I ε,κ ,j∈J ε,κ K ij -κ log(εκ -1 ) i∈I ε,κ µ i -κ -1 log(εκ) j∈J ε,κ ν j .
The above problem uses only the restricted parts K (Iε,κ,Jε,κ) , K (Iε,κ,J ε,κ ) , and K (I ε,κ ,Jε,κ) of the Gibbs kernel K for calculating the objective function Ψ ε,κ . Hence, a gradient descent scheme will also need only those rows/columns of K. This is in contrast to Sinkhorn algorithm which performs alternating updates of all rows and columns of K. In summary, SCREENKHORN consists of two steps: the first one is a screening pre-processing providing the active sets I ε,κ , J ε,κ . The second one consists in solving Equation ( 6) using a constrained L-BFGS-B [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF] for the stacked variable θ = (u Iε,κ , v Jε,κ ). Pseudocode of our proposed algorithm is shown in Algorithm 1. Note that in practice, we initialize the L-BFGS-B algorithm based on the output of a method, called RESTRICTED SINKHORN (see Algorithm 2 in the supplementary), which is a Sinkhorn-like algorithm applied to the active dual variables θ = (u Iε,κ , v Jε,κ ). While simple and efficient, the solution of this RESTRICTED SINKHORN algorithm does not satisfy the lower bound constraints of Problem (6) but provide a good candidate solution. Also note that L-BFGS-B handles box constraints on variables, but it becomes more efficient when these box bounds are carefully determined for problem (6). The following proposition (proof in supplementary material) expresses these bounds that are pre-calculated in the initialization step of SCREENKHORN. Proposition 1. Let (u sc , v sc ) be an optimal pair solution of problem (6) and K min = min i∈Iε,κ,j∈Jε,κ

K ij .
Then, one has

ε κ ∨ min i∈Iε,κ µ i ε(m -m b ) + max j∈Jε,κ νj nεκKmin m b ≤ e u sc i ≤ max i∈Iε,κ µ i mεK min , (7)
and εκ ∨ min j∈Jε,κ ν j ε(n -n b ) + κ max i∈Iε,κ µi mεKmin n b ≤ e v sc j ≤ max j∈Jε,κ ν j nεK min (8)
for all i ∈ I ε,κ and j ∈ J ε,κ .

Theoretical analysis and guarantees

This section is devoted to establishing theoretical guarantees for SCREENKHORN algorithm. We first define the screened marginals µ sc = B(u sc , v sc )1 m and ν sc = B(u sc , v sc ) 1 n . Our first theoretical result, Proposition 2, gives an upper bound of the screened marginal violations with respect to 1 -norm. Proposition 2. Let (u sc , v sc ) be an optimal pair solution of problem (6). Then one has

µ -µ sc 2 1 = O n b cκ + (n -n b ) C ∞ η + m b √ nmcµν K 3/2 min + m -m b √ nmKmin + log √ nm m b c 5/2 µν (9) Algorithm 1: SCREENKHORN(C, η, µ, ν, n b , m b )
Step 1: Screening pre-processing

1. ξ ← sort(µ r(K)), ζ ← sort(ν c(K)); //(decreasing order) 2. ε ← (ξ n b ζ m b ) 1/4 , κ ← ζ m b /ξ n b ; 3. I ε,κ ← {i = 1, . . . , n : µ i ≥ ε 2 κ -1 r i (K)}, J ε,κ ← {j = 1, . . . , m : ν j ≥ ε 2 κc j (K)}; 4. µ ← min i∈Iε,κ µ i , μ ← max i∈Iε,κ µ i , ν ← min j∈Jε,κ ν i , ν ← max j∈Jε,κ ν i ; 5. u ← log ε κ ∨ µ ε(m-m b )+ε∨ ν nεκK min m b , ū ← log μ mεKmin ; 6. v ← log εκ ∨ ν ε(n-n b )+ε∨ κ μ mεK min n b , v ← log ν nεKmin ; 7. θ ← stack(ū1 n b , v1 m b ), θ ← stack(u1 n b , v1 m b );
Step 2: L-BFGS-B solver on the screened variables

8. u (0) ← log(εκ -1 )1 n b , v (0) ← log(εκ)1 m b ; 9. û, v ← RESTRICTED SINKHORN(u (0) , v (0) ), θ (0) ← stack(û, v); 10. θ ← L-BFGS-B(θ (0) , θ, θ); 11. θ u ← (θ 1 , . . . , θ n b ) , θ v ← (θ n b +1 , . . . , θ n b +m b ) ; 12. u sc i ← (θ u ) i if i ∈ I ε,κ and u i ← log(εκ -1 ) if i ∈ I ε,κ ; 13. v sc j ← (θ v ) j if j ∈ J ε,κ and v j ← log(εκ) if j ∈ J ε,κ ; 14. return B(u sc , v sc ). and ν -ν sc 2 1 = O m b c 1 κ + (m -m b ) C ∞ η + n b √ nmcµν K 3/2 min + n -n b √ nmKmin + log √ nm n b c 5/2 µν , (10)
where c z = z -log z -1 for z > 0 and c µν = µ ∧ ν with µ = min i∈Iε,κ µ i and ν = min j∈Jε,κ ν j .

Proof of Proposition 2 is presented in supplementary material and it is based on first order optimality conditions for problem (6) and on a generalization of Pinsker inequality (see Lemma 2 in supplementary).

Our second theoretical result, Proposition 3, is an upper bound of the difference between objective values of SCREENKHORN and dual of Sinkhorn divergence (2). Proposition 3. Let (u sc , v sc ) be an optimal pair solution of problem (6) and (u , v ) is the pair solution of dual of Sinkhorn divergence (2). Then we have

Ψ ε,κ (u sc , v sc ) -Ψ(u , v ) = O R( µ -µ sc 1 + ν -ν sc 1 + ω κ ) .
where

R = C ∞ η + log (n∨m) 2 nmc 7/2 µν and ω κ = |1 -κ| µ sc 1 + |1 -κ -1 | ν sc 1 + |1 -κ| + |1 -κ -1 |.
Proof of Proposition 3 is exposed in the supplementary material. Comparing to some other analysis results of this quantity, see for instance Lemma 2 in [START_REF] Dvurechensky | Computational optimal transport: Complexity by accelerated gradient descent is better than by Sinkhorn's algorithm[END_REF] and Lemma 3.1 in [START_REF] Lin | On efficient optimal transport: An analysis of greedy and accelerated mirror descent algorithms[END_REF], our bound involves an additional term ω κ (with ω 1 = 0). To better characterize ω κ , a control of the 1 -norms of the screened marginals µ sc and ν sc are given in Lemma 2 in the supplementary material.

Numerical experiments

In this section, we present some numerical analyses of our SCREENKHORN algorithm and show how it behaves when integrated into some complex machine learning pipelines.

Setup

We have implemented our SCREENKHORN algorithm in Python and used the L-BFGS-B of Scipy.

Regarding the machine-learning based comparison, we have based our code on the ones of Python Optimal Transport toolbox (POT) [START_REF] Flamary | POT: Python optimal transport library[END_REF] and just replaced the sinkhorn function call with a screenkhorn one. We have considered the POT's default SINKHORN stopping criterion parameters and for SCREENKHORN, the L-BFGS-B algorithm is stopped when the largest component of the projected gradient is smaller than 10 -6 , when the number of iterations or the number of objective function evaluations reach 10 5 . For all applications, we have set η = 1 unless otherwise specified.

Analysing on toy problem

We compare SCREENKHORN to SINKHORN as implemented in POT toolbox1 on a synthetic example. The dataset we use consists of source samples generated from a bi-dimensional gaussian mixture and target samples following the same distribution but with different gaussian means. We consider an unsupervised domain adaptation using optimal transport with entropic regularization. Several settings are explored: different values of η, the regularization parameter, the allowed budget n b n = m b m ranging from 0.01 to 0.99, different values of n and m. We empirically measure marginal violations as the norms µµ sc 1 and νν sc 1 , running time expressed as TSINKHORN TSCREENKHORN and the relative divergence difference | C, P -C, P sc |/ C, P between SCREENKHORN and SINKHORN, where P = ∆(e u )K∆(e v ) and P sc = ∆(e u sc )K∆(e v sc ). Figure 2 summarizes the observed behaviors of both algorithms under these settings. We choose to only report results for n = m = 1000 as we get similar findings for other values of n and m.

SCREENKHORN provides good approximation of the marginals µ and ν for "high" values of the regularization parameter η (η > 1). The approximation quality diminishes for small η. As expected µµ sc 1 and νν sc 1 converge towards zero when increasing the budget of points. Remarkably marginal violations are almost negligible whatever the budget for high η. According to computation gain, SCREENKHORN is almost 2 times faster than SINKHORN at high decimation factor n/n b (low budget) while the reverse holds when n/n b gets close to 1. Computational benefit of SCREENKHORN also depends on η with appropriate values η ≤ 1. Finally except for η = 0.1 SCREENKHORN achieves a divergence C, P close to the one of Sinkhorn showing that our static screening test provides a reasonable approximation of the Sinkhorn divergence. As such, we believe that SCREENKHORN will be practically useful in cases where modest accuracy on the divergence is sufficient. This may be the case of a loss function for a gradient descent method (see next section).

Integrating SCREENKHORN into machine learning pipelines

Here, we analyse the impact of using SCREENKHORN instead of SINKHORN in a complex machine learning pipeline. Our two applications are a dimensionality reduction technique, denoted as Wasserstein Discriminant Analysis (WDA), based on Wasserstein distance approximated through Sinkhorn divergence [START_REF] Flamary | Wasserstein discriminant analysis[END_REF] and a domain-adaptation using optimal transport mapping [START_REF] Courty | Optimal transport for domain adaptation[END_REF], named OTDA.

WDA aims at finding a linear projection which minimize the ratio of distance between intra-class samples and distance inter-class samples, where the distance is understood in a Sinkhorn divergence sense. We have used a toy problem involving Gaussian classes with 2 discriminative features and 8 noisy features and the MNIST dataset. For the former problem, we aim at find the best two-dimensional linear subspace in a WDA sense whereas for MNIST, we look for a subspace of dimension 20 starting from the original 728 dimensions. Quality of the retrieved subspace are evaluated using classification task based on a 1-nearest neighbour approach.

Figure 3 presents the average gain (over 30 trials) in computational time we get as the number of examples evolve and for different decimation factors of the SCREENKHORN problem. Analysis of the quality of the subspace have been deported to the supplementary material (see Figure 6), but we can remark a small loss of performance of SCREENKHORN for the toy problem, while for MNIST, accuracies are equivalent regardless of the decimation factor. We can note that the minimal gains are respectively 2 and 4.5 for the toy and MNIST problem whereas the maximal gain for 4000 samples is slightly larger than an order of magnitude.

For the OT based domain adaptation problem, we have considered the OTDA with 1 2 ,1 group-lasso regularizer that helps in exploiting available labels in the source domain. The problem is solved using a majorization-minimization approach for handling the non-convexity of the problem. Hence, at each iteration, a SINKHORN/SCREENKHORN has to be computed and the number of iteration is sensitive to the regularizer strength. As a domain-adaptation problem, we have used a MNIST to USPS problem in which features have been computed from the first layers of a domain adversarial neural networks [START_REF] Ganin | Domain-adversarial training of neural networks[END_REF] before full convergence of the networks (so as to leave room for OT adaptation). Figure 4 reports the gain in running time for 2 different values of the group-lasso regularizer hyperparameter, while the curves of performances are reported in the supplementary material. We can note that for all the SCREENKHORN with different decimation factors, the gain in computation goes from a factor of 4 to 12, without any loss of the accuracy performance.

Conclusion

The paper introduces a novel efficient approximation of the Sinkhorn divergence based on a screening strategy. Screening some of the Sinkhorn dual variables has been made possible by defining a novel constrained dual problem and by carefully analyzing its optimality conditions. From the latter, we derived some sufficient conditions depending on the ground cost matrix, that some dual variables are smaller than a given threshold. Hence, we need just to solve a restricted dual Sinkhorn problem using an off-the-shelf L-BFGS-B algorithm. We also provide some theoretical guarantees of the quality of the approximation with respect to the number of variables that have been screened. Numerical experiments show the behaviour of our SCREENKHORN algorithm and computational time gain it can achieve when integrated in some complex machine learning pipelines.

7 Supplementary material

Proof of Lemma 1

Since the objective function Ψκ is convex with respect to (u, v), the set of optima of problem ( 3) is non empty. Introducing two dual variables λ ∈ R n + and β ∈ R m + for each constraint, the Lagrangian of problem (3) reads as

L (u, v, λ, β) = ε κ λ, 1n + εκ β, 1m + 1 n B(u, v)1m -κu, µ - v κ , ν -λ, e u -β, e v
First order conditions then yield that the Lagrangian multiplicators solutions λ * and β * satisfy

∇uL (u * , v * , λ * , β * ) = e u * (Ke v * -λ * ) -κµ = 0n,
and ∇vL (u

* , v * , λ * , β * ) = e v * (K e u * -β) - ν κ = 0m
which leads to

λ * = Ke v * -κµ e u * and β * = K e u * -ν κe v *
For all i = 1, . . . , n we have that e u * i ≥ ε κ . Further, the condition on the dual variable λ * i > 0 ensures that e u * i = ε κ and hence i ∈ I ε,κ . We have that λ * i > 0 is equivalent to e u * i ri(K)e v * j > κµi which is satisfied when ε 2 ri(K) > κµi. In a symmetric way we can prove the same statement for e v * j .

Proof of Proposition 1

We prove only the first statement ( 7) and similarly we can prove the second one ( 8). For all i ∈ Iε,κ, we have e u sc i > ε κ or e u sc i = ε κ . In one hand, if e u sc i > ε κ then according to the optimality conditions λ sc i = 0, which implies e u sc i m j=1 Kije v sc j = κµi. In another hand, we have

e u sc i min i,j Kij m j=1 e v sc j ≤ e u sc i m j=1 Kije v sc j = κµi.
We further observe that m j=1 e v sc j = j∈Jε,κ e v sc j + j∈J ε,κ e v sc j ≥ εκ|Jε,κ| + εκ|J ε,κ | = εκm. Then max i∈Iε,κ e u sc i ≤ maxi∈I ε,κ µi mεKmin .

Analogously, one can obtain for all j ∈ Jε,κ

max j∈Jε,κ e v sc j ≤ maxj∈J ε,κ νj nεKmin . (11)
Now, since Kij ≤ 1, we have

e u sc i m j=1 e v sc j ≥ e u sc i m j=1 Kije v sc j = κµi.
Using (11), we get

m j=1 e v sc j = j∈Jε,κ e v sc j + j∈J ε,κ e v sc j ≤ εκ|J ε,κ | + maxj∈J ε,κ νj nεKmin |Jε,κ|.
Therefore,

min i∈Iε,κ e u sc i ≥ ε κ ∨ κ minI ε,κ µi εκ(m -m b ) + max j∈Jε,κ ν j nεK min m b .

Proof of Proposition 2

We define the distance function :

R+ × R+ → [0, ∞] by (a, b) = b -a + a log( a b ).
While is not a metric, it is easy to see that is not nonnegative and satisfies (a, b) = 0 iff a = b. The violations are computed through the following function:

d (γ, β) = n i=1 (γi, βi), for γ, β ∈ R n + .
Note that if γ, β are two vectors of positive entries, d (γ, β) will return some measurement on how far they are from each other. The next Lemma is from Abid and Gower (2018) (see Lemma 7 herein).

Lemma 2. For any γ, β ∈ R n + , the following generalized Pinsker inequality holds γ -

β 1 ≤ 7( γ 1 ∧ β 1)d (γ, β).
The optimality conditions for (u sc , v sc ) entails

µ sc i = e u sc i m j=1 Kije v sc j , if i ∈ Iε,κ, ε κ m j=1 Kije v sc j , if i ∈ I ε,κ = κµi, if i ∈ Iε,κ, ε κ m j=1 Kije v sc j , if i ∈ I ε,κ , (12)
and

ν sc j = e v sc j n i=1 Kije u sc i , if j ∈ Jε,κ, εκ n i=1 Kije u sc i , if j ∈ J ε,κ = ν j κ , if j ∈ Jε,κ, εκ n i=1 Kije u sc i , if j ∈ J ε,κ . (13)
By ( 12), we have

d (µ, µ sc ) = n i=1 µ sc i -µi + µi log µi µ sc i = i∈Iε,κ (κ -1)µi -µi log(κ) + i∈I ε,κ ε κ m j=1 Kije v sc j -µi + µi log µi ε κ m j=1 Kije v sc j = i∈Iε,κ (κ -log(κ) -1)µi + i∈I ε,κ ε κ m j=1 Kije v sc j -µi + µi log µi ε κ m j=1 Kije v sc j .
Now by ( 8), we have in one hand

i∈I ε,κ ε κ m j=1 Kije v sc j = i∈I ε,κ ε κ j∈Jε,κ Kije v sc j + εκ j∈J ε,κ Kij ≤ i∈I ε,κ ε κ m b max i,j Kij maxj∈J ε,κ νj nεKmin + (m -m b )εκ max i,j Kij ≤ (n -n b ) m b maxj νj nκKmin + (m -m b )ε 2 .
On the other hand, we get

ε κ m j=1 Kije v sc j = ε κ j∈Jε,κ Kije v sc j + εκ j∈J ε,κ Kij ≥ m b Kmin mε 2 Kmin minj∈J ε,κ νj κ((n -n b )mε 2 Kmin + mε 2 Kmin + n b κ maxi∈I ε,κ µi) + ε 2 (m -m b )Kmin ≥ mm b ε 2 (Kmin) 2 minj∈J ε,κ νj κ((n -n b )mε 2 Kmin + mε 2 Kmin + n b κ maxi∈I ε,κ µi) + ε 2 (m -m b )Kmin ≥ mm b ε 2 K 2 min minj∈J ε,κ νj κ((n -n b )mε 2 Kmin + mε 2 Kmin + n b κ maxi∈I ε,κ µi) .
Then

1 ε κ m j=1 Kije v sc j ≤ κ((n -n b )mε 2 Kmin + mε 2 Kmin + n b κ maxi∈I ε,κ µi) mm b ε 2 K 2 min minj∈J ε,κ νj ≤ κ(n -n b + 1) m b Kmin minj∈J ε,κ νj + n b κ 2 maxi∈I ε,κ µi mm b ε 2 K 2 min minj∈J ε,κ νj . It entails i∈I ε,κ ε κ m j=1 Kije v sc j -µi + µi log µi ε κ m j=1 Kije v sc j ≤ (n -n b ) m b nκKmin + (m -m b )ε 2 -min i µi + max i µi log κ(n -n b + 1) maxi µi m b Kmin minj∈J ε,κ νj + n b κ 2 (maxi µi) 2 mm b ε 2 K 2 min minj∈J ε,κ νj . Therefore d (µ, µ sc ) ≤ n b cκ max i µi + (n -n b ) m b maxj νj nκKmin + (m -m b )ε 2 -min i µi + max i µi log κ(n -n b + 1) maxi µi m b Kmin minj∈J ε,κ νj + n b κ 2 (maxi µi) 2 mm b ε 2 K 2 min minj∈J ε,κ νj .
Finally, by Lemma 2 we obtain

µ -µ sc 2 1 ≤n b cκ max i µi + 7(n -n b ) m b maxj νj nκKmin + (m -m b )ε 2 -min i µi + max i µi log κ(n -n b + 1) maxi µi m b Kmin minj∈J ε,κ νj + n b κ 2 (maxi µi) 2 mm b ε 2 K 2 min minj∈J ε,κ νj .
Following the same lines as above, we also have

ν -ν sc 2 1 ≤m b c 1 κ max i µi + 7(m -m b ) n b κ maxi µi mKmin + (n -n b )ε 2 -min j νj + max j νj log (m -m b + 1) maxj νj n b κKmin mini∈I ε,κ µi + m b (maxj νj) 2 nn b ε 2 κ 2 K 2 min mini∈I ε,κ µi .
To get the closed forms ( 9) and ( 10), we used the following facts: Remark 1. We have log(1/K r min ) = r C ∞/η, for every r ∈ N. Using (5), we further derive:

ε = O((mnK 2 min ) -1/4 ), κ = O( m/(ncµν Kmin)), κ -1 = O( n/(mKmincµν ), (κ/ε) 2 = O(m 3/2 / √ nKmin(cµν ) 3/2 ), and (εκ) -2 = O(n 3/2 / √ mKminc 3/2 µν ).
1.

Next, we bound the two terms u sc -u ∞ and v sc -v ∞. If r ∈ I ε,κ , then we have

|(u sc )r -u r | = log m j=1 Krje v j m j=1 κµr mε ( ) ≤ log max 1≤i≤m Krje v j κµr mε ≤ max 1≤j≤m (v j -log( κµr mε ) ≤ v -log( κµr mε ) ∞ ≤ v -v sc ∞ + log( mε 2 cµν ).
where the inequality ( ) comes from the fact that n j=1 a j n j=1 b j ≤ max 1≤j≤n a j b j , ∀aj, bj > 0. Now, if r ∈ Iε,κ, we get

|u sc r -u r | = log κ m j=1 Krje v j m j=1 Krje (v sc ) j ≤ log m j=1 Krje v j m j=1 Krje (v sc ) j ( ) ≤ v sc -v ∞. If s ∈ J ε,κ then |v sc s -v s | = log(εκ) -log( νs n i=1 Kise u i ) ≤ log max 1≤i≤n Kise u i νs nκε ( ) ≤ max 1≤i≤n (u i -log( νs nκε ) ≤ u -log( νs nκε ) ∞ ≤ u -u sc ∞ + log( nε 2 cµν ).
If s ∈ Jε,κ then

|v sc s -v s | = log κ m i=1 Krie u i m i=1 Krie (u sc ) i ≤ log κ m i=1 Krie v j m i=1 Krie (u sc ) i ( ) ≤ u sc -u ∞.
Therefore, we obtain the followoing bound:

max{ u -u sc ∞, v -v sc ∞} ≤ max u ∞ + u sc ∞ + log( nε 2 cµν ), v ∞ + v sc ∞ + log( mε 2 cµν ) ≤ 2 u ∞ + v ∞ + u sc ∞ + v sc ∞ + log (n ∨ m)ε 2 cµν . (15)
Now, Lemma 3.2 in [START_REF] Lin | On efficient optimal transport: An analysis of greedy and accelerated mirror descent algorithms[END_REF] provides an upper bound for the ∞ of the optimal solution pair (u , v ) of problem (2) as follows: u ∞ ≤ A and v ∞ ≤ A, where

A = C ∞ η + log n ∨ m c 2 µν . (16)
Plugging ( 15) and ( 16) in ( 14), we obtain

Ψε,κ(u sc , v sc ) -Ψ(u , v ) ≤ 2 A + u sc ∞ + v sc ∞ + log (n ∨ m)ε 2 cµν ) µ sc -µ 1 + ν sc -ν 1 + |1 -κ| u sc ∞ µ sc 1 + µ sc -µ 1 (17) + |1 -κ -1 | v sc ∞ ν sc 1 + ν sc -ν 1 . By Proposition 1, we have u sc ∞ ≤ log ε κ ∨ 1 mεKmin and v sc ∞ ≤ log εκ ∨ 1 nεKmin
and hence by Remark 1,

u sc ∞ = O log(n 1/4 /(mKmin) 3/4 c 1/4 µν ) and u sc ∞ = O log(m 1/4 /(nKmin) 3/4 c 1/4 µν ) .
Acknowledging that log(1/K 2 min ) = 2 C ∞/η, we have

A + u sc ∞ + v sc ∞ + log (n ∨ m)ε 2 cµν ) = O C ∞ η + log (n ∨ m) 2 nmc 7/2 µν . Letting Ωκ := |1 -κ| u sc ∞ µ sc 1 + µ sc -µ 1 + |1 -κ -1 | v sc ∞ ν sc 1 + ν sc -ν 1 . We have that Ωκ = O C ∞ η + log 1 (nm) 3/4 c 1/2 µν |1 -κ|( µ sc 1 + µ sc -µ 1 + |1 -κ -1 |( ν sc 1 + ν sc -ν 1) = O C ∞ η + log (n ∨ m) 2 nmc 7/2 µν |1 -κ| µ sc 1 + |1 -κ -1 | ν sc 1 + |1 -κ| + |1 -κ -1 | .
Hence, we arrive at

Ψε,κ(u sc , v sc ) -Ψ(u , v ) = O R( µ -µ sc 1 + ν -ν sc 1 + ωκ) .
To more characterize ωκ, the following lemma expresses an upper bound with respect to 1-norm of µ sc and ν sc .

Lemma 3. Let (u sc , v sc ) be an optimal solution of problem (6). Then one has

µ sc 1 = O n b √ m nKmincµν + (n -n b ) m b √ nmcµν K 3/2 min + m -m b √ nmKmin , (18)
and

ν sc 1 = O m b √ n mKmincµν + (m -m b ) n b √ nmcµν K 3/2 min + n -n b √ nmKmin . ( 19
)
Proof. Using inequality (8), we obtain

µ sc 1 = i∈Iε,κ µ sc i + i∈I ε,κ µ sc i (12) = κ µ sc Iε,κ 1 + ε κ i∈I j∈Jε,κ Kije v sc j + εκ j∈J ε,κ Kij (8) ≤ κ µ sc Iε,κ 1 + (n -n b ) m b maxj∈J ε,κ νj nκKmin + (m -m b )ε 2 .
Using Remark 1, we get the desired closed form in ( 18). Similarly, we can prove the same statement for ν sc 1.

Additional experimental results

Experimental setup. All computations have been run on each single core of an Intel Xeon E5-2630 processor clocked at 2.4 GHz in a Linux machine with 144 Gb of memory.

On the use of a constrained L-BFGS-B solver. It is worth to note that standard Sinkhorn's alternating minimization cannot be applied for the constrained screened dual problem ( 6). This appears more clearly while writing its optimality conditions (see Equations ( 12) and ( 13) ). We resort to a L-BFGS-B algorithm to solve the constrained convex optimization problem on the screened variables ( 6), but any other efficient solver (e.g., proximal based method or Newton method) could be used. The choice of the starting point for the L-BFGS-B algorithm is given by the solution of the RESTRICTED SINKHORN method (see Algorithm 2), which is a Sinkhorn-like algorithm applied to the active dual variables. While simple and efficient the solution of this RESTRICTED SINKHORN algorithm does not satisfy the lower bound constraints of Problem (6). We further note that, as for the SINKHORN algorithm, our SCREENKHORN algorithm can be accelerated using a GPU implementation2 of the L-BFGS-B algorithm [START_REF] Fei | Parallel L-BFGS-B algorithm on GPU[END_REF].

Algorithm 2: RESTRICTED SINKHORN

1. set: fu = εκ c(K Iε,κ,J ε,κ ), fv = εκ -1 r(K I ε,κ ,Jε,κ ); 2. for t = 1, 2, 3 do f (t) v ← K Iε,κ,Jε,κ u + fv ; v (t) ← ν Jε,κ κf (t) v ; f (t) u ← K Iε,κ,Jε,κ v + fu ; u (t) ← κµ Iε,κ f (t) u ; u ← u (t) , v ← v (t) ; 3. return (u (t) , v (t) )
Comparison with other solvers. We have considered experiments with GREENKHORN algorithm [START_REF] Altschuler | Near-linear time approximation algorithms for optimal transport via Sinkhorn iteration[END_REF] but the implementation in POT library and our Python version of Matlab Altschuler's Greenkhorn code3 were not competitive with SINKHORN. Hence, for both versions, SCREENKHORN is more competitive than GREENKHORN. The computation time gain reaches an order of 30 when comparing our method with GREENKHORN while SCREENKHORN is almost 2 times faster than SINKHORN. We note a slight loss of performance for the toy problem, whereas for MNIST, all approaches yield the same performance. (bottom-right). Gain in computation time. This toy problem is a problem in which classes are overlapping and distance between samples are rather limited. According to our analysis, this may be a situation in which SCREENKHORN may result in smaller computational gain. We can remark that with respect to the accuracy SCREENKHORN with decimation factors up to 10 are competitive with SINKHORN, although a slight loss of performance. Regarding computational time, for this example, small decimation factors does not result in gain. However for above 5-factor decimation, the gain goes from 2 to 10 depending on the number of samples. The accuracy panels also report the performance of a 1-NN when noadaptation is performed. We remark that the strenght of the class-based regularization has influence on the performance of SCREENKHORN given a decimation factor. For small value on the left, SCREENKHORN slightly performs better than SINKHORN, while for large value, some decimation factors leads to loss of performances. Regarding, running time, we can note that SINKHORN is far less efficient than SCREENKHORN with an order of magnitude for intermediate number of samples. chosen among the set {0.1, 1, 5, 10}. We can note that in this situation, there is not loss of accuracy while our SCREENKHORN is still about an order of magnitude more efficient than Sinkhorn.

  we denote the complementary set of L by L = {1, . . . , n}\L. We also denote |L| the cardinality of L. Given a vector w ∈ R n , we denote w L = (w i1 , . . . , w i k ) ∈ R k and its complementary w L ∈ R n-k . The notation is similar for matrices; given another subset of indices S = {j 1 , . . . , j l } ⊆ {1, . . . , m} with j 1 < • • • < j l , and a matrix T ∈ R n×m , we use T (L,S) , to denote the submatrix of T , namely the rows and columns of T (L,S) are indexed by L and S respectively. When applied to matrices and vectors, and (Hadamard product and division) and exponential notations refer to elementwise operators. Given two real numbers a and b, we write a ∨ b = max(a, b) and a ∧ b = min(a, b).

Figure 1 :

 1 Figure 1: Plots of (e u , e v ) with (u , v ) is the pair solution of dual of Sinkhorn divergence (2) and the thresholds α u , α v .

  This guarantees that |I ε,κ | = n b and |J ε,κ | = m b by construction. In addition, when (n b , m b ) tends to the full number budget of points (n, m), the objective in problem (3) converges to the objective of dual of Sinkhorn divergence (2).

Figure 2 :

 2 Figure 2: Empirical evaluation of SCREENKHORN vs SINKHORN for normalized cost matrix i.e. C ∞ = 1. (most-lefts): marginal violations in relation with the budget of points on n and m . (center-right) ratio of computation times TSINKHORN TSCREENKHORN and, (right) relative divergence variation. The results are averaged over 30 trials.

Figure 3 :Figure 4 :

 34 Figure 3: Wasserstein Discriminant Analysis : running time gain for (left) a toy dataset and (right) MNIST as a function of the number of examples and the data decimation factor in SCREENKHORN.

  Figure 5: TGREENKHORN TSCREENKHORN : Running time gain for the toy problem (see Section 5.2) as a function of the data decimation factor in SCREENKHORN, for different settings of the regularization parameter η.

Figure 6 :

 6 Figure 6: Accuracy of a 1-nearest-neighbour after WDA for the (left) toy problem and, (right) MNIST).We note a slight loss of performance for the toy problem, whereas for MNIST, all approaches yield the same performance.

Figure 7 :Figure 8 :

 78 Figure7: (top-left) Accuracy and (bottom-right) computational time gain on the toy dataset for η = 0.1 and 1-nearest-neighbour. (bottom) accuracy and gain but for a 5-nearest-neighbour. We can note that a slight loss of performances occur for larger training set sizes especially for 1-nearestneighbour. Computational gains increase with the dataset size and are on average of the order of magnitude.

Figure 9 :

 9 Figure9: OT Domain adaptation MNIST to USPS : (top) Accuracy and (bottom) running time of SINKHORN and SCREENKHORN for hyperparameter of the p,1 regularizer (left) λ = 1 and (right) λ = 10. Note that this value impacts the ground cost of each Sinkhorn problem involved in the iterative algorithm. The accuracy panels also report the performance of a 1-NN when noadaptation is performed. We remark that the strenght of the class-based regularization has influence on the performance of SCREENKHORN given a decimation factor. For small value on the left, SCREENKHORN slightly performs better than SINKHORN, while for large value, some decimation factors leads to loss of performances. Regarding, running time, we can note that SINKHORN is far less efficient than SCREENKHORN with an order of magnitude for intermediate number of samples.

Figure 10 :

 10 Figure10: OT Domain adaptation MNIST to USPS : (left) Accuracy and (right) running time of SINKHORN and SCREENKHORN for the best performing (on average of 10 trials) hyperparameter p,1 chosen among the set {0.1, 1, 5, 10}. We can note that in this situation, there is not loss of accuracy while our SCREENKHORN is still about an order of magnitude more efficient than Sinkhorn.

https://pot.readthedocs.io/en/stable/index.html

https://github.com/nepluno/lbfgsb-gpu

https://github.com/JasonAltschuler/OptimalTransportNIPS17

Acknowledgments

This work was supported by grants from the Normandie Projet GRR-DAISI, European funding FEDER DAISI and OATMIL ANR-17-CE23-0012 Project of the French National Research Agency (ANR).

Proof of Proposition 3

We first define K a rearrangement of K with respect to the active sets Iε,κ and Jε,κas follows:

Setting .

µ = (µ Iε,κ , µ I ε,κ

) , .

ν = (ν Jε,κ , ν J ε,κ

) and for each vectors u ∈ R n and v ∈ R m we set .

u = (u Iε,κ , u I ε,κ

) and .

) . We then have

v, and

v, where

Let us consider the convex function

v sc ) 1n . Gradient inequality of any convex function g at point xo reads as g(xo) ≥ g(x) + ∇g(x), xo -x , for all x ∈ dom(g). Applying the latter fact to the above function at point (u , v ) we obtain 1n, B( .

Then,

Applying Holder's inequality gives

where, in the last inequality, we use the facts that κ .

Then Ψε,κ(u sc , v sc ) -Ψ(u , v ) ≤ u sc -u ∞ + |1 -κ| u sc ∞ µ sc -µ 1