Sparsistency for inverse optimal transport
Authors: Francisco Andrade
Gabriel Peyré
Clarice Poon
Venue: ICLR 2024 poster (2024)
Keywords: Optimal transport
sparsity
sparsistency
metric learning
Primary area: metric learning, kernel learning, and sparse coding

Abstract: Optimal Transport is a useful metric to compare probability distributions and to compute a pairing given a ground cost. Its entropic regularization  variant (eOT) is crucial to have fast algorithms and reflect fuzzy/noisy matchings. This work focuses on Inverse Optimal Transport (iOT), the problem of inferring the ground cost from samples drawn from a coupling that solves an eOT problem. It is a relevant problem that can be used to infer unobserved/missing links, and to obtain meaningful information about the structure of the ground cost yielding the pairing. On one side, iOT benefits from convexity, but on the other side, being ill-posed, it requires regularization to handle the sampling noise. This work presents an in-depth theoretical study of the $\ell_1$ regularization to model for instance Euclidean costs with sparse interactions between features.  Specifically, we derive a sufficient condition for the robust recovery of the sparsity of the ground cost that can be seen as a far reaching generalization of the Lasso’s celebrated ``Irrepresentability Condition’’. To provide additional insight into this condition (consequently on the types of recoverable costs) we work out in detail the Gaussian case. Surprisingly, varying the entropic regularizer provides evidence that the Gaussian iOT interpolates between a graphical Lasso and a classical Lasso, thereby establishing a connection between iOT and graph estimation, an important problem in ML.

Published as a conference paper at ICLR 2024
SPARSISTENCY FOR INVERSE OPTIMAL TRANSPORT
Francisco Andrade & Gabriel Peyré
Clarice Poon
ENS Paris
University of Warwick clarice waweeepooonwwwaaww .k
(f (franci sco. andrade, gabriel .eterriee.pereeeeee Bens. fr
ABSTRACT
Optimal Transport is a useful metric to compare probability distributions and to compute pairing given ground cost. Its entropic regularization variant (eOT) crucial to have fast algorithms and reflect fuzzy/noisy matchings. This work focuses on Inverse Optimal Transport (iOT), the problem of inferring the ground cost from samples drawn from a coupling that solves an eOT problem. It is relevant problem that can be used to infer unobserved/ missing links, and obtain meaningful information about the structure of the ground cost yielding the pairing. On one side, iOT benefits from convexity, but on the other side, being ill ill-posed, it requires regularization to handle the sampling noise. This work presents an in-depth theoretical study of the é regularization to model for instance Euclidean costs with sparse interactions between features. Specifically, we derive a sufficient condition for the robust recovery of the sparsity of the ground cost that can be celebrated "Irrepresentability Condition". To provide additional insight into this condition, work out in detail the Gaussian case. We show that as the entropic penalty varies, the iOT problem interpolates between a graphical Lasso and classical Lasso, thereby establishing connection between iOT and graph estimation, an important problem in ML.
far reaching generalization of the Lasso's
INTRODUCTION
Optimal transport has emerged key theoretical and numerical ingredient in machine learning for performing learning over probability distributions. It enables the comparison of probability distributions in "geometrically faithful" manner by lifting a ground cost (or "metric" in a loose sense) between pairs of points to a distance between probability distributions, metrizing the convergence in law. However, the success this OT approach hypothesis that the ground cost is adapted to the problem under study. This necessitates the exploration of ground metric learning. However, it is exceptionally challenging due to its priori highly non-convex nature when framed complications in its mathematical analysis. As we illustrate this theoretical article, these problems become tractable numerically and theoretically- if one assumes access to samples from the OT coupling (i.e., having access to some partial matching driven by the ground cost), Admittedly, this is restrictive setup, but it arises in practice (refer subsequent sections for illustrative applications) and can also be construed as a step in a more sophisticated learning pipeline. The purpose of this paper is to propose some theoretical understanding of the possibility of stably learning a ground cost from partial matching observations.
inherently tied the
optimization problem, thereby inheriting
1.1 PREVIOUS WORKS
Entropic Optimal Transport. OT has been instrumental defining and studying various procedures at the core of many ML pipelines, such as bag- of-feetures matching Rubner et al. (2000), distances in NLP Kusner et al. (2015), generative modeling Arjovsky et al. (2017), flow evolution for sampling Bortoli al. [2021), , and even single-cell trajectory inference SeneernceeSchiee Schiebinger lttal, 2019). We refer to the monographs Santambrogio (2015) for detailed accounts on the theory OT, and Peyré al. (2019) for its computational aspects. Of primary importance to work, entropic regularization of OT is the workhorse of many ML applications. It enables a fast and highly parallelizable estimation of the OT coupling using the Sinkhorn algorithm Sinkhorn (1964).

---

Published as a conference paper at ICLR 2024
More importantly, it defines smooth distance that incorporates the understanding that matching procedures should be modeled as a noisy process (i.e., should not be assumed to be 1:1). These advantages were first introduced in ML by the seminal paper of Cuturi Cuturi (2013) and this approach finds its roots in Schrödinger's work in statistical physics Lénnard 2012). The role of noise in matching (with applications in economics) and its relation entropic OT were advanced in series of papers by Galichon and collaborators Galichon & Salanié (2010); Dupuy & Calichon (2014; ; Galichon & Salanić (2022); ; see the book Galichon (2018). These works are key inspirations
for the present paper, which aims at providing more theoretical understanding in the case of inverse OT (as detailed next).
Metric Learning. The estimation of some metrics from pairwise interactions (either positive or negative) falls into the classical field of metric learning in ML, and we refer to the monograph Bellet et al. (013)) for more details. In contrast to the inverse OT (iOT) problem considered in this paper, classical metric learning is more straightforward. as no global matching between sets of points is involved. This allows the metric to be directly optimized, while the iOT problem necessitates some form of bilevel optimization. Similarly to our approach, since the state space is typically continuous, it is necessary to restrict the class of distances to render the problem tractable. The common option, which we also adopt in our paper to exemplify our findings, is to consider the class of Mahalanobis distances. These distances generalize the Euclidean distance and are equivalent to computing a vectorial embedding of the data points. See, for instance, Xing et al. 2002 Weinberger al.
2006 Davis & Dhillon 2008
OT Ground Metric Learning. The problem of estimating the ground cost driving OT in a super- vised manner was first addressed (utuuri & Avis (014). Unlike methods that have access to pairs of samples, the ground metric learning problem requires pairs of probability distributions and then evolves into a classical metric learning problem, but within the OT space. The class of ground metrics can be constrained, for example, by utilizing Mahalanobis Wang & Guibas 2012); Xu et al. 2018); Kerdoneuff et al. (002)) geodesic distances Heitz et al. (021)), to devise more efficient learning schemes. The BrudyZen et al. (0014) conducts ground metric learning and matrix factorization simultaneously, finding applications
NLP Huang al. 20166) It is noteworthy that ground metric learning can also be linked to generative models through adversarial train- ing[Gnnvvay et al. (2018) and to robust learning & Cuturi (2020 by maximizing the cost to render the OT distance as discriminative as possible.
Inverse Optimal Transport. The inverse optimal transport problem (IOT) can be viewed as a specific instance of ground metric learning, where one aims to infer the ground cost from partial observations of the (typically entropically regularized) optimal transport coupling. This problem was first formulated and examined by Dupuy and Galichon Dupuy & Galichon 2014) over discrete space (also see Galichon & Salanié (002) for a more detailed analysis), making the fundamental remark that the maximum likelihood estimator amounts to solving convex problem. The mathematical properties of the iOT problem for discrete space (i.c., direct computation of the cost between all pairs of points) are explored in depth in Chiu et al. (2022), studying uniqueness (up to trivial ambiguities) and stability to pointwise noise. Note that our theoretical study differs fundamentally we focus continuous state spaces. This "continuous" setup assumes access only to a set of couples, corresponding to matches (or links) presumed to be drawn from an OT coupling. In this scenario, the iOT is typically an ill ill-posed problem, and (Dupuy et al. (2019); Carlier et al. (2023) propose regularizing the maximum likelihood estimator with either a low- rank (using nuclear norm penalty) or sparse prior (using an e Lasso-type penalty). In our work, we concretely focus on the sparse case, but our theoretical treatment of the iOT could be extended to general structured convex regularization, along the lines of Vaiter et al. 2015). While not the focus of our paper, it is noteworthy that these works also propose efficient large-scale, non-smooth proximal solvers to optimize the penalized maximum likelihood functional, and we refer Ma et al. (2020) for an efficient solver without inner loop calls to Sinkhorn's algorithm. This approach was further refined in Stuart & Wolfram (2020), deriving it from a Bayesian interpretation, enabling the use of MCMC methods to sample the posterior instead of optimizing pointwise estimate (as we consider here). They also propose parameterizing cost functions as geodesic distances graphs (while consider only linear models to maintain convexity). An important application of iOT to ML, explored Li al. (2019), is to perform link prediction by solving new OT problems once the cost has been estimated from the observed couplings. Another category
2

---

Published as a conference paper at ICLR 2024
of ML application of iOT is representation learning (learning embeddings of data into, e.g., Euclidean space) from pairwise interactions, as demonstrated by Shi et al. 20233 which recasts contrastive learning specific instance iOT.
Inverse problems and model selection. The iOT problem is formally a bilevel problem, the ob- servation model necessitates solving an OT problem as an inner-level program program(Colson Colson et al. (2005)
-we refer to Eisenberger et al. (2022) for recent numerical treatment of bilevel programming with entropic OT. The iOT can thus be conceptualized as an "inverse optimization" problem ZZhng Liu 1996); Ahuja Orlin (2001), but with a particularly favorable structure, allowing it to be recast as a convex optimization. This provides the foundation for a rigorous mathematical analysis of performance, as we propose in this paper. The essence of our contributions a theoretical examination of the recoverability of the OT cost from noisy observations, and particularly, the robustness to noise of the sparse support of the cost (for instance, viewed as symmetric matrix for Mahalanobis norms). There exists rich tradition of similar studies in the fields of inverse problem regularization and model selection in statistics. The most prominent examples are the sparsistency theory of the Lasso Tibshirani 19966 (least square regularized by (t), which culminated in the theory of compressed sensing Candès et al.] (2006). These theoretical results are predicated on so- co-called "irrepresentability condition" Zhao & Yu (2006), which ensures the stability of the support. While our analysis grounded in similar concepts (in particular, we identify the corresponding irrepresentability condition for the iOT), the iOT inverse problem is fundamentally distinct due to the differing observation model (it corresponds to a sampling process rather than the observation of vector) and the estimation process necessitates solving linear program of potentially infinite dimension (in the limit of large number of samples). This mandates a novel proof strategy, which forms the core of our mathematical analysis.
1.2 CONTRIBUTIONS
This paper proposes the first mathematical analysis of the performance regularized iOT esti- mation, focusing on the special case of sparse é regularization (the e- -iOT method), We begin by deriving the customary "irrepresentability condition" of the iOT problem, rigorously proving that it is well-defined. This condition interweaves the properties of the Hessian of the maximum likelihood functional with the sparse support of the sought-after cost. The main contribution of this paper is Theorem 5 which leverages this abstract irrepresentability condition to ensure sparsistency the iio method. This relates to the robust estimation of the cost and its sup- port in some linear model, assuming the number of samples is large enough. Specifically, demonstrate a sample complexity of n 1/22 Our subsequent sets of contributions are centered on the case of matching between samples of Gaussian distributions. Herein, we illustrate in Lemma] how to compute the irrepresentability condition in closed form. This facilitates the examination of how the parameters of the problem, particularly regularization strength and the covariance of the distributions, influence the success and stability of iOT. We further explore the limiting cases of small and large entropic regularization, revealing in Proposition @ and Proposition that iOT interpolates between the graphical lasso (to estimate the graph structure of the precision matrix) and a classical lasso. This sheds light on the connection between iOT and graph estimation procedures. Simple synthetic numerical explorations in Section further provide intuition about how and the geometry of the graph associated with sparse cost impact sparsistency. As a minor numerical contribution, present in Appendix F large-scale el -iOT solver, implemented in JAX and distributed as open- aper-source software.
2 INVERSE OPTIMAL TRANSPORT
The forward problem Given probability distributions a (2)), B ( and cost function c: R, the entropic optimal transport problem seeks to compute a coupling density
Sink(c, arggmax(c, π) E -KL(naa @ () where (c) clx, y)dT(x,)) (a,f)
(1)
where KL(T() & fllggdd/ddddd f dn and V (a, B) is the space of all probability measures (. x 2) with marginals f.
3

---

Published as a conference paper at ICLR 2024
The inverse problem The inverse optimal transport problem seeks recover the cost function given an approximation An of the probability coupling f Sink(c, A typical setting is an empirical probability coupling n = - oun,y) where (x, Vo) iid ft. See Section 2.2
The loss functton The iOT problem has been proposed and studied in series of papers, see Section 1.2 The approach is typically to consider some linear parameterization of the cost CA(X, y) by some parameter A R2. The key observation of Dupuy et al. 2019 that the negative log-likelihood of i at parameter value A is given by
E (A, _ - (CA, ft + (4(A) where (*(A) & sup (CA, π) KL(Tlå @ P),
2
(
and å, β are the marginals of ft. For ease of notation, unless stated otherwise, we write W W So, computing the maximum likelihood estimator A for the cost corresponds minimizing the convex lloss function' A (A, R), which, by regarding Wj as convex conjugate, can seen as an instance of Venchel-Yourgy loss, family of losses proposed in Blondel al. 2020) We write the parameterization as
Φ A - =A= Σ A,C,, where (2 F=1
A relevant example are quadratic loss functions, so that for 2 c c 3 c given A e Rdixd, () x Ay. In this case, s did2 and for k= = (f, j) E (d1] x (d2), Cx y) x;y].
4-10T To handle the presence of noisy data (typically coming from the sampling process), various regularization approaches have been proposed. In this work, we focus on the use of e-regulaarization Carlier et al. 2023 recover sparse parametrizations:
argmin& where S 2 ||||I + ((A,
$
(iOT-((
Kantorovich formulation Note that ((A)) is defined via concave optimization problem and by Fenchel duality, one can show that (OT has the following equivalent Kantorovich formulation Carlier et al. (0033:
argmin H (A.f.g), where ((A.f.)) ((A.f.) A1IA and A.f.g
(Hoo)
(A. -- ( y 2. g()) da(x)dp(y).
Based on this formulation, various algorithms have been proposed, including alternating mini- mization with proximal updates Carlier et al. (2023 Section details Fettalls new large scale solver that we use for the numerical simulations.
2.1 INVARIANCES AND ASSUMPTIONS
Assumptton 1. We first assume that I and y are compact.
Note that y has the translation invariance property that for any constant function u, f(f + u. u, A) (f. A), so to remove this invariance, throughout, we restrict the optimization of [
to the set
y a ())
(a) x
122())) 8(y)dB()=00
(2)
Next, we make some assumptions on the cost to remove invariances in the iOT problem. Assumption 2 (Assumption on the cost).
(i) Exxx)- y) z Id is invertible.
iii) CCxx, y & 1 for a almost every and (-almost every y.
In this work, we restrict linear parameterizations, although the same loss function can also be applied learn with nonlinear parameterization, e.g. via neural networks Ma 2020).
4

---

Published as a conference paper at ICLR 2024
(iii) for all k, ('Cg(x, for f.a.e. and y)dp(y) for a-a.e. X.
Under these assumptions, it can be shown that iOT has unique solution (see remark after Propo- sition Assumption (i) is to ensure that сд = @A is uniquely determined by Assumption 2 (ii) is without loss of generality, since assume that are compactly supported, so this holds up to a rescaling of the space. Assumption 2 (iii) is to handle the invariances pointed out in Carlier et al. 2023) and Ma al. (020)):
I (A.f.g)) I [ f'g)) CA gf c (f' eg').
As observed in Carlier et al. 2023 any - ( Letting C₂(x, y) C2(x, y)
be adjusted to fit this assumption: can define
where (xxxx) [x(x and yx(y) / (x,((x)))))()))
4 = Σx AcC we have [If $ -8--2. (g + (f DA. So optimization with the parametrization
equivalent to optimization with
NB: For the quadratic cost (x(x,y) = x,yj for k= conditio (iii) corresponds to recentering the data points, and taking x x- xxlx(x and y y- ydß(y). Condition (ii) holds if [|x|] v || for a-a.e. and f-a.e. y. Condition corresponds to invertibility
EplyyT].
2.2 THE FINITE SAMPLE PROBLEM
In practice, we do not observe f but data-points (x) ld To recover A, plug into i the empirical measure An estimator
A for i 1,..., n1 where 2 = Sink(c E). 1 n f 6x,,, and consider the
argmin ||A f).
(iOT-(1 (x),
As section 2.1 to account for the invariance when solving fOT - (xt we again centre the cost parameterization such that for all Σ Cg and for all j. Σラ- also that ( ()) can be formulated entirely in finite dimensions. See Appendix Bor details.
Cx vj) 0. Note
THE CERTIFICATE FOR SPARSISTENCY
In this section, consider the problem ( with full data ft and present a sufficient condition for support recovery, that we term non-degeneracy of the certificate. For simplicity of notation, throughout this section denote W & Wg. Under non-degen non-degeneracy of the certificate we obtain support recovery stated in Theorem 3 known result whose proof can be found e.g. Lee et al. (2015). This condition can be seen as a generalization of the celebrated Lasso's Irrepresentability condition (see e.g. (Hattie et al. )2015) Lasso corresponds to having quadratic loss instead of (A, thus V2 (AA) in the definition below reduces to a matrix. In what follows, we denote uj (uj)ier and U1,J & (Ui,))icI.je) the restriction operators. Definition The certificate with respect to A and support Ai *
(V2 (" sign(A))
(C)
We say that it is non-degenerate if (= Iloo <.
The next proposition, whose proof can be found Appendix C.T shows that the function W(A) twice differentiable, thus ensuring that well defined.
Proposition A A- ((A) is twice differentiable, strictly convex, with gradient and Hessian
W V W(A) = θπ
(x.
aA
where n, A is the unique solution to 7 with cost CA AAA
We remark that strict convexity of W implies that any solution of (1T must be unique, and by I T-convergence, solutions A to 1iT ( converge to Å as À - 0. The next theorem,
well-known result (see (Lee et al. 2015 Theorem 3.4) for a proof), shows that support recovery can be characterized via 2;
5

---

Published as a conference paper at ICLR 2024
Theorem 3. Leti Sink(eg.E). non-deggeneaee,, then for all A sufficientll small, the solution A to ( is sparsistent with || || -Al (A) and Supp (A*) Supp(A).
3.1 INTUITION BEHIND THE CERTIFICATE
Implication of the non-degenyracy condition The non-degeneracy condition widely studied for the Lasso problem Hastie et al. 2015 Just as for the Lasso, the closer (z²)) is to the more
> then one cannot expect to recover this coefficients when there is noise. Note also that in the Lasso, the pre-certificate formula in C roughly the correlation between coefficient inside and outside the support I. our case, the loss is more complex and this correlation is measured according to the hessian of W, which integrates the curvature of the loss. This curvature formula however involved, so to gain intuition, we perform detailed analysis in the and E- co for the Gaussian setting where it becomes much simpler (see Section 5.1
unstable is this coefficient, and (z)))
Link to optimality conditions The certificate 2 can be seen the limit optimality condition for the optimization problem ( (?) as 0: by the first order optimality condition to ( A is solution if and only if 2 è [VL(A*) A where the subdifferential for the ( norm has the explicit form |||A| = {z: 2z 1, Vi Supp((), = sign(A,)) It follows z can be seen as certificate for the support of A since Supp(A*) = {i: |z = 1/. study the support behavior of A for small A, it is therefore interesting to consider the limit of x as À - 0. Its limit is precisely the subdifferential element with the minimal norm and coincides with ( under the nondegeneracy condition:
Proposition 4. Let 2
VI(AA) where A solves 10 Then,
lim -² 2 A argmin ((z, (V" 2)p 2E Alli]
(MNC)
1-0
Moreover, 2 is non-degenerate, then min A
2
4 SAMPLE COMPLEXITY BOUNDS
Our main contribution shows that a certificate for sparsistency under sampling noise: Theorem 5. Let Sink(cq.E). Suppose that the certificate 2; non-degenerate. Let 0. Then, for all sufficiently small regularization parameters and sufficiently many number of samples n,
A and max(exp(C||]])) /E) (og(11/))== V log(2s)
S nn,
for some constant C 0, with probability at least - 6, the minimizer An to
? is sparsistent with
À with Supp(An) Supp(A) and A - 5 A+ exp(C |A]] /1)
Main idea behind Theorem 5 We know from Theorem 3 that there is some An 0 such that for all & λa, the solution to % has the same support as Å. show that the finite sample problem also recovers the support of Á when is sufficiently large, we fix A E (0, Ag] and consider the setting where the observations are iid samples from the coupling measure ft We will derive convergence bounds for the primal and dual solutions the number of samples increases. Let (Aor fow. (oo) minimise
K.o
Denote Foo = (foo (x))))))) Go (go0 (y)) and
where pro(x, y) exp (x, y) foolx) + +go()))
Poo n (x)
Let An minimize (0T-f1(A) (f)). Then, there exists probability matrix Pn E R and vectors Fn, G, Rº" such that P, + exp [ (0,, F2 @ Gn)). In fact, (An, Fn, G₂) minimize the finite dimensional dual problem K given in the appendix. Now consider the certificates
- 0" (poc 8 8- and Zn A
(P, -P))
2o0
1
6

---

Published as a conference paper at ICLR 2024
Note that Zoo and Zn both depend on A; the superscript was dropped since is fixed. Moreover, Zoo is precisely ZA from Propostion 4 By exploiting strong convexity properties of fn one can show the following sample complexity bound the convergence of 2n Zoo (the proof can be found the appendix):
Proposition 6. Let n Z max (log(1 ( ..gg((2)) for some > O. For some constant C 0, with probability at least 1 - 6, 1o 2n 2 exp(C) All, /e))og(1/8)A and
AA, - 2 (Fnn - n
- G exp(C|| ( /
From Proposition for all A Åo for some Åo sufficiently small, the magnitude of 2f outside the support of less than one. Moreover, the convergence result in Proposition that, for n sufficiently large, the magnitude of în outside the support of Å is also less than one. Hence, since the set (i; (zn); = +1 determines the support of Ap, we have sparsistency.
implies
5 GAUSSIAN DISTRIBUTIONS
To get further insight about the sparsistency property of iOT, we consider the special case where the source and target distributions are Gaussians, and the cost parametrization Ca y) = Ay. To this end, we first derive closed form expressions for the Hessian d², ( = V (A)). Given .V (ma Σa) and = V (mp, Ep), it is known (see Bojilov & Galichon 2016)) that the coupling
a Σ
density is also Gaussian of the form M can be written an optimization problem over the cross-oovariance x BBojilov Galichon (2016).
(m
for some Σ Rd xdz In this case, W Σp
Σ
W(A) ΣER ddl
sup (A, Σ) logdet(Ip - Σ >>=
(4)
with the optimal solution being precisely the cross-covariann of the optimal coupling x. In Bojilov & Galichon 20166), the authors provide an explicit formula for the minimizer Σ to and, consequently, for ((A):
g2
Σ == A=AAAN(AA -- --tT where
A
EB+
AtI A
(5)
4
By differentiating the first order condition for W, that is At derives the expression for the Hessian in terms of 5
' (EB - Σ ΣΣZ ")'') Galichon
E
A V = -, = =
(A 1.T & 'A
(6)
where such that Tvec(A) vec(AT). This formula does not hold when rectangular or rank deficient, since At is not differentiable. In the following, we derive, via the implicit function theorem, a general formula for 0, that agrees with that of Galichon the square invertible case.
Lemma Denoting Σ as in 5), one has
= -
(A e + eAT)"
This formula given in Lemma provides an explicit expression for the certificate
C
5.1 LIMIT CASES FOR LARGE AND SMALL
E
This section explores the behaviour of the certificate in the large/ limits: Proposition reveals that the large e limit coincides with the classical Lasso while 9 9 rrevelss that the small epsilon limit (for symmetric x and Σ, = Σg = Id) coincides with the Graphical Lasso. In the following results, we denote the functional in fioT ( with parameters A and & by F (A).
Proposition 8 (E - oo). Let À be invertible and let Ωt Sink(c be the observed coupling between = (m,, Ex) and ß = N (mg, Σz). Then,
lim 2 (E, x En) ann (EE & Ealunn) sign(A)r. 000
(8)
7

---

Published as a conference paper at ICLR 2024
Moreover, for Ao O, given any (ejnuenceeee))) and A argmin e, (A) with lim 00, any cluster point of [Aj) is in
1 argmin Ao Aoll + 11= 1/12/2 @ 2/ (A- - A
(9)
AeRded
Interpretation E -SE-OO, the KL term forces Å to be close to the independent coupling a @ f and in the limit, iOT is simply Lasso problem and the limit in 8 is precisely the Lasso certificate Hastie et al. (2015). Here, the cross covariance of ft satisfies e A (e) so for large E, sparsity in A indicates the independence in the coupling between a and B.
53 in the appendix),
Proposition 9 (s 0). Let Ả be symmetric positive-definite and let á Sink(cg, E) be the observed coupling between a = N (ma, Id) and B = N (mg, Id). Then,
lim 4. = ( 1 s ; Jen 7 (a,n) sign(A)r- 0-
(10)
Let Ao > 0. Then, optimizing over symmetric positive semi- refindet matrices, given any sequence f and Aj argmin Azo Azo Fcooe, (A) with j-m é 1 = 0, any cluster point of (Aj) is in
argmin λa All logdet(A)) + (A. 7-1 0=0
(11)
2
Interpretation In contrast, E - O, the KL. term disappears and the coupling f becomes de- pendent. Naturally, the limit problem is the graphical lasso, typically used infer conditional independence in graphs (but where covariates can be highly dependent). Note also that the limit (10) precisely the graphical Lasso PrtificateHHastie et al. 2015 Here, (Remark (55 in the appendix) one show that for (x, y) π, the conditional covariance of x conditional on (and also vice versa) EA- 0 (e²2) Sparsity in A can therefore be viewed as information on conditional independence.
5.2 NUMERICAL ILLUSTRATIONS
In order to gain some insight into the impact o and the covariance structure on the efficiency of iOT, present numerical computations of certificates here. We fix the covariances of the input measures as Σn Σ3 Id,, similar results are obtained with different covariance as long as they are not rank-deficient. We consider that the support of the sought-after cost matrix A 8ld,, diag(GI,)) Ge "** is defined as shifted Laplacian matrix of some graph adjacency matrix G, for graph of size n = 80 (similar conclusions hold for larger graphs). We set the shift 0 to be 10% of the largest eigenvalue of the Laplacian, ensuring that C is symmetric and definite. This setup corresponds to graphs defining positive interactions at vertices and negative interactions along edges. For small adopting the graphical lasso interpretation (as exposed in Section 5.1 and interpreting C as a precision matrix, this setup corresponds (for instance, for a planar graph) to imposing spatially smoothly varying covariance C-1, Figure illussratees how the value of the certificates ž1,j evolves depending on the indexes (i, j) for three types of graphs (circular, planar, and Erdós-Rént construction, Z;, = and Zi.j = -1 for (i. f) connected by the graph. For 2 to be non degeeegeeer,,,, it is required that | < 1, moves away from j on the graph. For the circular and planar graphs, the horizontal axis represents the geodesic distance (i. demonstrating how the certificates become well-behaved the distance increases. The planar graph displays envelope curves showing the range of values of z for a fixed value of dgeod j), while this is a single-valued curve for the circular graph due to periodicity. For the Erdös-Rénn Mends-Rényieee graph, to better account for the randomness of the certificates along the graph edges, we display the histogram of the distribution of i for dgeod(i, (i. j)= = 2 (which represents the most critical set of edges, as they are the most likely to be large). All these examples show the same behavior, namely, that increasing improves the behavior of the certificates (which is in line with the stability analysis of Sectionn5.1 since (8) implies that for large e, the certificate is trivially non -eggeeeeeee whenever Σạ Σn diagonal), and that pairs of vertices (i, j) connected by a small distance dgeod (i, j) are the most likely to be degenerate. This suggests that they will be inaccurately estimated by iOT for small
with probability of edges equal to 0.1), for several values of E. By
8

---

Published as a conference paper at ICLR 2024
65
68-
-10
Circular
Planar
Erd Rényi
Figure 1: Display of the certificate values zı,j for three types of graphs, for varying Left, middle:
plotted as a function of the geodesic distance dgeod(i,) (i, j) on the x-axis. Right: histogram of 21, for (i distance dgeod(i, j) =
Figure 2 displays the recovery performances P'-iOT for the circular graph shown on the left of Figure (similar results are obtained for the other types of graph topologies). These numerical simulations are obtained using the large-scale iOT solver, which we detail in Appendix The performance is represented using the number of inaccurately estimated coordinates in the esti- mated cost matrix C, so a score of 0 means a perfect estimation of the support (sparsistency is achieved). For 0.1, sparsistency cannot be achieved, aligning with the fact that the certificate is degenerate as depicted in the previous figure. In sharp contrast, for larger &, sparsistency can be attained as soon the number of samples N is sufficiently large, which also aligns with our theory of sparsistency of (--iOT since the certificate is guaranteed to be non non-degenerate in the large setting.
c = 0.1
E=1
== 10
V-165
wrongly estim. pos.
6x2266
22
29
50
063
C.M
0.59
016
0.56
0.55
0.26
0.00
Figure Recovery performance (number of wrongly estimated position) (1-iOT as function A for three different values of e.
Figure 3 in the appendix displays the certificate values in the case of a non-symmetric planar graph. The graph is obtained by deleting all edges on planar graph with < We plot the certificate values as function of geodesic distances (i, (i,
of the symmetrized graph. The middle plot shows the certificate values on > (where the actual edges are constrained to be and nondegeneracy requires values smaller than in absolute value for dgeod (i. /) > 2). The right plot shows the certificate values on & where there are no edges and for nondegeneracy. one expects values smaller than in absolute value for dgeod j) 1. Observe that here, the certificate degenerate for small values of (i (i, j) when = 0, meaning that the problem unstable at the "ghost" symmetric edges. As co, the certificate becomes non-degenerate.
9

---

Published as a conference paper at ICLR 2024
CONCLUSION
In this paper, we have proposed the first theoretical analysis of the recovery performance of -iOT. Much this analysis can be extended to more general convex regularizers, such as the nuclear norm to promote low-rank Euclidean costs, for instance. Our analysis and numerical exploration support the conclusion that iOT becomes ill -ol-po and fails maintain sparsity for overly small E. When approached from the perspective of graph estimations, unstable indices occur at smaller geodesic distances, highlighting the geometric regularity of iOT along the graph geometry.
ACKNOWLEDGEMENTS
The work of G. Peyré was supported by the European Research Council (ERC project NORIA) and the French government under management of Agence Nationale de la Recherche as part of the "Investissements d' d'avenir" program, reference ANR19-P3IA-0001 (PRAIRIE 3IA Institute).
REFERENCES
Ravindra K Ahuja and James B Orlin. Inverse optimization. Operations research, 49(5):771-783, 2001.
Jason Altschuler, Jonathan Niles- and Philippe Rigollet. Near- linear time approximation algo- rithms for optimal transport via sinkhorn iteration. Advances in neural information processing systems, 30, 2017.
Martin Arjovsky, Soumith Chintala, and Léon Bottou. Wasserstein generative adversarial networks. In International conference on machine learning, pp. 214-223. PMLR, 2017.
Aurélien Bellet, Amaury Habrard, and Marc Sebban. A survey on metric learning for feature vectors and structured data. arXiv preprint arXiv:1306.6700 2013.
Mathieu Blondel, André FT Martins, and Vlad Niculae. Learning with fenchel-young losses. The Journal of Machine Learning Research, 21(1):1314-1382, 2020.
Raicho Bojilov and Alfred Galichon. Matching in flosed-form: equilibrium, identification, and comparative statics. Economic Theory, 61:587-609, 2016.
Jonathan Borwein and Adrian Lewis. Convex Analysis. Springer, 2006.
Emmanuel J Candès, Justin Romberg, and Terence Tao. Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information. IEEE Transactions on information theory, 52(2):489 509, 2006.
Guillaume Carlier, Arnaud Dupuy, Alfred Galichon, and Yifei Sun. Sista: learning optimal transport costs under sparsity constraints. Communications on Pure and Applied Mathematics, 76(9): 1659-1677, 2023.
Wei-Ting Chiu, Pei Wang, and Patrick Shafto. Discrete probabilistic inverse optimal transport. In International Conference on Machine Learning, pp.
PMLR, 2022.
Benoît Colson, Patrice Marcotte, and Gilles Savard. Bilevel programming: survey. 4or, 3:87-107, 2005.
Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. In Adu in Neural Information Processing Systems, pp. 2292-2300, 2013.
Marco Cuturi and David Avis. Ground metric learning. The Journal of Machine Learning Research, 15(1):533-564, 2014.
Marco Cuturi and Gabriel Peyré. A smoothed dual approach for variational wasserstein problems. SIAM Journal on Imaging Sciences, 9(1):320-343 2016.
10

---

Published as a conference paper at ICLR 2024
Jason V Davis and Inderjit S Dhillon. Structured metric learning for high dimensional problems. In Proceedings of the 14th ACM SIGKDD international conference on Knowledge discovery and data mining, pp. 195-208, 2008.
Valentin De Bortoli, James Thorttn, Jeremy Heng, and Arnaud Doucet. Diffusion schrödinger bridge with applications to score-based generative modeling. Advances in Neural Information Processing Systems, 14:16695-77709 2021.
Klaus Deimling. Nonlinear vonlinearuunctioall analysis. Courier Corporation, 2010.
Armaud Dupuy and Alfred Galichon. Personality traits and the marriage market. Journal of Political Economy, 122(6):1271-1113 2014.
Arnaud Dupuy, Alfred Galichon, and Yifei Sun. Estimating matching affinity matrices under low-rank constraints. Information and Inference: Journal of the IMA, 8(4):677-689, 2019.
Marvin Eisenberger, Aysim Toker, Laura Leixe Florian Bernard, and Daniel Cremers. unified framework for implicit sinkhorn differentiation, In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 509-518, 2022.
Alfred Galichon. Optimal transport methods economics. Princeton University Press, 2018.
Alfred Galichon and Bernard Salanié. Matching with trade-offs: Revealed preferences over com- peting characteristics. Preprint hal-00473173, 2010.
Alfred Galichon and Bernard Salanié. Cupid's invisible hand: Social surplus and identification in matching models. The Review of Economic Studies, 89(5):2600-2629, 2022.
Aude Genevay, Gabriel Peyré, and Marco Cuturi. Learning generative models with sinkhorn divergences. In International Conference on Artificial Intelligence and Statistics, pp. 1608-1617.
PMLR, 2018.
Aude Genevay, Lénaicc Chizat, Francis Bach, Marco Cuturi, and Gabriel Peyré. Sample complexity of Sinkhorn divergences. In The 22nd international conference on artificial intelligence and statistics, pp. 1574-1583, PMLR, 2019.
Trevor Hastie, Robert Tibshirani, and Martin Wainwright. Statistical learning with sparsity: the lasso and generalizations. CRC press, 2015.
Matthieu Heitz, Nicolas Bonneel, David Cocurjolly, Marco Cuturi, and Gabriel Peyré. Ground metric learning on graphs. Journal of Mathematical Imaging and Vision, 63:89-107, 2021.
Jean-Baptiste Hiriart- and Claude Lemaréchal. Conjugacy in convex analysis. Convex Analysis and Minimization Algorithms II: Advanced Theory and Bundle Methods, pp. 35-90, 1993.
Gao Huang, Chuan Guo, Matt J Kusner, Sun, Fei Sha, and Kilian Q Weinberger. Supervised word mover's distance. Advances in neural information processing systems, 29, 2016.
Tanguy Kerdoncuff, Rémi Emonet, and Mare Sebban. Metric Learning in Optimal Transport for
Domain Adaptation. In International Joint Conference on Artificial Intelligence, Kyoto, Japan, January 2021. URL URL httphtttp//haala..
archi ves ves www..ives-.varrrrs..f.ror.e ft/mjrrresrr//ujmmmmee -200:mmmm226111
Matt Kusner, Yu Sun, Nicholas Kolkin, and Kilian Q Weinberger. From word embeddings to document distances. In Proc. of the 32nd Intern. Conf. on Machine Learning, pp. 957-966, 2015.
Jason D Lee, Yuekai Sun, and Jonathan E Taylor. On model selection consistency of regularized m-estimators. Electronic Journal of Statistics, 9:608-642, 2015.
Christian Léonard. From the Schrödinger problem to the Kantoroantoo problem. Journal of Functional Analysis, (4889)4)88791199999 2012.
Ruilin Li, Xiaojing Ye, Haomin Zhou, and Hongyuan Zha. Learning to match via inverse optimal transport. Journal of machine learning research, 20, 2019.
11

---

Published as a conference paper at ICLR 2024
Shaojun Ma, Haodong Sun, Xiaojing Ye, Hongyuan Zha, and Haomin Zhou. Learning cost functions for optimal transport. arXiv preprint rrXiv:2002.09650 2020.
Gonzalo Mena and Jonathan Niles-Wer Weed. Statistical bounds for entropic optimal transport: sample complexity and the central limit theorem. Advances in neural information processing systems, 32, 2019.
Marcel Nutz. Introduction to entropic optimal transport. Lecture notes, Columbia University, 2021.
François-Pierre Paty and Marco Cuturi. Regularized optimal transport is ground cost adversarial. In International Conference on Machine Learning, pp. 7532-7542. PMLR, 2020.
Gabriel Peyré, Marco Cuturi, et al. Computational optimal transport: With applications to data science. Foundations and Trends@ in Machine Learning, 11(6-6):55--07,,22 2019.
Clarice Poon and Gabriel Peyré. Smooth bilevel programming for sparse regularization. Advances in Neural Information Processing Systems, 34:1543-1s 1555,, 2021.
Philippe Rigollet and Austin Stromme. On the sample complexity of entropic optimal transport. arXiv preprint arXiv:2206. 2tti::2220....7.2.. 13472, 2022.
Yossi Rubner, Carlo Tomasi, and Leonidas J. Guibas. The earth mover's distance as metric for image retrieval. International Journal of Computer Vision, 40(2)99--21, November 2000.
Filippo Santambrogio. Optimal transport for applied mathematicians. Birkäuser, NY, 55(58-63):94, 2015.
Geoffrey Schiebinger, Jian Shu, Marcin Tabaka, Brian Cleary, Vidya Subramanian, Aryeh Solomon, Joshua Gould, Siyan Liu, Stacie Lin, Peter Berube, et al. gene expression identifies developmental trajectories reprogramming. Cell, 176(4):928 943, 2019.
analysis of single-cell
Liangliang Shi, Gu Zhang, Haoyu Zhen, Jintao Fan, and Junchi Yan. Understanding and generaliz- ing contrastive learning from the inverse optimal transport perspective. 2023.
Richard Sinkhorn. relationship between arbitrary positive matrices and doubly stochastic matrices. Ann. Math. Statist., 35:876-879, 1964.
Andrew M Stuart and Marie-Theres, Wolfram. Inverse optimal transport. SIAM Journal on Applied Mathematics, 80(1):599-611 2020.
Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society Series B: Statistical Methodology, 58(1):667-888, 1996.
Joel A Tropp et al. An introduction to matrix concentration inequalities. Foundations and Trends@ in Machine Learning, 8(1-2):1-230, 2015.
Samuel Vaiter, Mohammad Golbabaee, Jalal Fadili, and Gabriel Peyré. Model selection with low complexity priors. Information and Inference: Journal of the IMA, 4(3):230-287, 2015.
Fan Wang and Leonidas J Guibas. Supervised earth mover's distance learning and its computer vision applications. In European Conference on Computer Vision, pp. 442-455. Springer, 2012.
Kilian Q Weinberger, John Blitzer, and Lawrence K Saul. Distance metric learning for large margin nearest neighbor classification. In Advances in neural information processing systems, pp. 1473-1480, 2006.
Eric Xing, Michael Jordan, Stuart Russell, and Andrew Ng. Distance metric learning with applica- tion to clustering with side-information, Advances in neural information processing systems, 15, 2002.
Jie Xu, Lei Luo, Cheng Deng, and Heng Huang. Multi-level metric learning via smoothed wasser- stein distance. In IJCAI, pp. 2919-29229 2018.
12

---

Published as a conference paper at ICLR 2024
Gloria Zen, Elisa Ricci, and Nicu Sebe. Simultaneous ground metric learning and matrix factoriza- tion with earth mover's distance. In 2014 22nd International Conference on Pattern Recognition, pp. 3690-3695, IEEE, 2014.
Jianzhong Zhang and Zhenhong Liu. Calculating some inverse linear programming problems. Journal of Computational and Applied Mathematics, 72(2):261-273) 1996.
Peng Zhao and Bin Yu. On model selection consistency of lasso. The Journal of Machine Learning Research, 7:2541-2563, 2006.
Figure 3: Display of certificate values for non-symmetric planar graph, for varying with edges only for j>} Middle/ /Right plots of the certificate values as a function the geodesic distance dgeodlt, j) of the symmetrized graph. The middle plot show the values when restricted to j.. The right plot shows the values restricted to & (where there are no edges).
A INTERPRETATIONS OF THE LOSS FUNCTION
As mentioned in Section the loss LAA, f) can be recovered via maximum likelihood (ML) estimation (Dpuyy al. 2019 Proposition 1) and be regarded as an instance of a family of losses called Fenchel- renchel-Young losses. For the sake of completeness, this section explains this in more detail.
A.1
MAXIMUM LIKELIHOOD INTERPRETATION
The map A Si, where Sink(.,)) is defined in Section 2 set of measures which absolutely continuous with respect to By the standard duality result (see Nutz (2021) )
can be seen as parameterizing
1) The density of Sink(CA,E) is given by d'iin(cc.
exp fA ga), where fA and gA solve the dual problem, iee,, sup J da I g expica f+ d(a @ () + +1;
2) The values of the primal and dual problems agree and equal to fA da ga dp.
Combining these two facts we obtain that
dSink(cA.E) CA 84
log d(a ())
sup (C - ne
KL(zia
W*(A),
where W3 is defined in as in Section Finally, taking expectation with respect to A yields
Ex.yy log dSink(cA.E)
Exx.y-- +Exx,n (Wa()) -c (4(A) L(A,))
d(a @ ()
thus establishing the connection with MI. estimation.
Bilevel interpretation
Note also that
argmin ((A, = argmin - (CA. -
2
argmin (CA. nA)
KI
_K s
2
13

---

Published as a conference paper at ICLR 2024
where A argmin (CA, πτ) KL(T) B). Since T, have the same marginals,
(c,. - nイ> = (cΛ + fa SA. TA
(log(T A/l @ B)), a).
We therefore have
argmin L (A arg (L
where
RA = argmin - B).
A
A
xe
A.2 FENCHEL FENCHEL-YOUNNG Loss INTERPRETATION
is argued in Blondel al. 2020 that, associated a prediction rule of the form
ý(0) argmax (8, 2(p),
pedom(Q)
there is natural loss function that the authors term Fenchel-Your Young loss and that is given by
() A Q (0) (2(p) (8.
With this in mind, let ((x) = KL e p) Lz (z) where (a)),,)))) (+) is the indicator function of (a, (), and note that
L ft) = La -
Since, for the purposes of minimization with respect to A, the term Q(+) is irrelevant, we see that finding a minimizer (T) amounts to finding a minimizer of the I regularized Fenchel-Young loss associated with KL.- @ (3) (a)
B THE FINITE SAMPLE PROBLEM
As mention in Section 2.2
we do not observe the full coupling á = Sink(cg,E), but only the the empirical measure n = Σ ōx,,; with (x,, y)) itd at showow here that the problem ( can be formulated entirely in finite dimensions follows. Let
H(PIQ) Σ Pl.goooPii
Note that n has marginals â, Σ 8x, and b, E We can interpret ft, as the matrix P, Idnxx and the "noisy" primal problem can be equivalently written as
min sup АLАD + (,,,, P Pn) ACR' PER
E H(P)
s.f.
P1=-1
1 and p = 11,
(9P)
2
where we write H(P) H(P) 1) and 4 Ar where C₁ = (C2 ma** Note that the dinte-dimensional problem has the same invariances as [IOT E (ft) so, we
will take Cg be centred that for all [,(Co), =0 and for all j. (C,)), The finite
sample Kantorovich formulation is
inf (,,AFF where (,(A,F IA(AF + A|AA|| and A.F.G.P
(Kr)
(0, A). (Pn) [F 2n Σ exp 2 (F: (F+ () (0,
I.(A,
and we restrict the optimization of H over Sr, ((A.F. &G 0
C PROOFS FOR SECTION
C.1 PROOF OF PROPOSITION DIFFERENTIABILITY OF W)
For the strict convexity of ((A) see Lemma 3 in Dupuy Galichon (001). The gradient formula can also be found in Dupuy & Galichon 2014) and trivially follows from the envelope theorem
14

---

Published as a conference paper at ICLR 2024
because the optimization problem W(A) has unique solution. We will only give proof of the Hessian formula.
The formula for the Hessian follows from the formula for the gradient provided we show that the density A is continuously differentiable with respect to A. The fact that we can swap the order of the operator @' and partial differentiation follows from the conditions for differentiability under the integral sign which holds since the measures are compactly supported. Without loss of generality, let 1. Then, the optimizer in ((A) is of the form Santambrogio (2015):
xA(x,y) = exp [A(x) + va(y) + cc(xx,)))
where u-(+) and va(+) satisfy
ux(x) = log
exp pa(y) +raly)+c((x,1)) apty), a-a.s.
va(y) = -log exp Ua(x) + ca(x, y)) da(x)) B-a.s.
It is known (see e.g. Nutz (021) that the functions (-(+) and v-(-) inherit the modulus of con- tinuity of the cost (in this case the map (x,y) - y))) and, hence, since we are assuming the measures to be compactly supported, it follows that U4 E (a) and VA E (f). Moreover, (uA, P') solve these two equations then, for any constant c, the pair (4 + C, V c) also solution and, hence, to eliminate this ambiguity we consider solutions in H L (a) (B)), where 12(00) = (p); f ddp()) 0}.
To show that W twice differentiable, since (W OD aT is sufficient show that A 44 and A - VΛ differentiable. To this end, we will apply the Implicit Function Theorem (in Banach spaces) to the map
F: R (c) (B)
u log for exp u(y) + cx(x, dBly) P lo1o exp ((x)) + cca(x, y) da(x))
P
A
since we have F(UA, VA = O. F(uA, UA, A), is the linear map defined by
The partial derivative of F at (AA, denoted by
(r+ PAl. [ (xlxx))()da()))) (12)
where
exp va(y)) ca(x))) and y)(x,
exp [a(x)) +
p)(x,))=
fø exp + dp(y)
Sx exp un(y) + ca(x, df(x)
Note that, since F(UA, PA, A) = PA qa(x,)) xx(x, Moreover, since A has marginals and B, follows that
S f(xxdd((xx++++++)) [ xx 2f d(a y) + Jyy gly) apty).
(13)
=f Jxxe (F(r) ((0++ al y) d(a (xx
This shows that Fu,A UA A) invertible - the last line 13 is zero if and only if gg= and since g (B) it follows that and f 0.
To conclude the proof we need to show that (Ou, A)) a bounded operator (see e.g. Deimling (2110) for the statement of IFT in Banach spaces) and to show this it is enough show that, for some constant C,
(Qu. Un A))f.. Cll
15

---

first show that the vector z coincides with the solution dual Published as a conference paper at ICLR 2024
This follows from 13 and the fact that there exists constant C such that xA(x, for all and (see e.g. Nutz
2021 for the existence of C). In fact, from g E (p),,, we obtain f(f(x) +
g(y))) d(a @ (x. II( g)l? and, hence, 13 implies that
(If.g) (Q.. F(uA, VA A)) (f) > cil(f
and from Cauchy-Schwarza appled to the left-hand-siiee we obtain
vA. ( g)ll.
thus concluding the proof.
C.2 CONNECTION OF THE CERTIFICATE
(C) WITH OPTIMIZATION PROBLEM 1OT (0
As mentioned in Section 3.1 the vector 2 provides insight into support recovery since Supp(4*) s {i; z = =11 In this section, we provide the proof to Proposition 4 which shows that as Å con- verges to 0, d converges to the solution of a quadratic optimization problem (MNC that term the minimal norm certificate.
Note that the connection with can now be established by noting that if the minimal norm certificate is non-degenerate then the inequality constraints (which correspond to the complement of the support of Ā) in [MC] can be dropped since they are inactive; in this case (MNC) reduces to quadratic optimization problem with only equality constraints and whose solution can be seen to be (C) which will thus be non-degenerate as well. The converse clear. So, under non- degeneracy, can be seen as the limit optimality vector and determines the support of A when A is small.
To prove Proposition problem of [OT-2 (7)).
Proposition 10. Let W be the convex conjugate W. Problem HOT admits a dual given by
argmin W (Éxy A2) subject to
2| < 1,
(14)
where 40 f. Moreover, pair of primal-dual solutions solutions Ah,2², ) is related by
a =
and Y' Al
(15)
Proof. Observe that can write ((A,
-(A,a) = W(A)
Ja (o ()x, di W(A) ( A)r A) W(A) - (A,E2y)
The Fenchel Duality Theorem (see egg. Borwein & Lewis 2006 ) yields a dual of (0T-(1 given by
argmin W W' (Exy- w) (A)) (1) (w).
To conclude the proof just note that the Fenchel conjugate of A is the indicator of the set
{v: || & A) and make change of variable z è 1/Aw to obtain 14 The relationship between any primal-dual pair in Fenchel Duality can also be found Borwein Lewis 2006).
Proof of Proposition 4 that
We begin by noting that W* is class C² in a neighborhood of Exy and
V (
(16)
To see this, note that Proposition together with the assumption on á implies that Σxy= VA
16

---

Published as a conference paper at ICLR 2024
Moreover, since ((A) is twice continuously differentiable and strictly convex (see Proposition 2), follows (see egg. Corollary 4.2.9 in Hiriartt- Urruty & Lemaréchal (1993)) that W (-) is C² and strictly convex in neighborhood of Exy and that 16) holds.
Now observe that, since VW (Lxy) vw (wW.(WW(A) A, we can rewrite d [|]]] as
o| | = argmin - VW (Exy)) subject to 2| $
(17)
Observe that since zamdda are uniformly bounded vectors due to the constraint set in there convergent subsequence converging to some z". We later deduce that all limit points are the same and hence, the full sequence z converges to z". Let An be such that ,, -0 = 2", and let 2 be any element in @|| 0||A|| We have that
<- 1-2.Vir VW (xy)) & (- ",,"=="''''" (xx))) & _ (W" (Exy-Anzzz) w w(Exy))
| nn
(w W An ,n - w(Er)))
where the first inequality is the optimality of 2", the second inequality is the gradient inequality of convex functions and the last inequality follows from the optimality of 2 Taking the limit as A, - 0 we obtain that
2 VW
ve'VV"'r)) (L,))) = = (-22, v (Exy)),
showing that 2" E0 Alp We now finish the proof by showing that
x v" (2", (Ee)22))
و v2
(18)
Since (*)-) is Cº in neighborhood of Σxy, Taylor's theorem ensures that there exists remainder function R(x) with lim,- R(x) = such that
An An, v (Ex))=n)) + R(Ê.y + A Anz2n) An VW (xxr)> 2
( An w (Exy))
1
W (2xy 1 Anz w (Exy)
An
An
An (20, w" (Ex R(Exy* + An )
=(-2), (
2
An (2
<(-z^> (xy)) "
v² ( + RlExy + ,, 11
2
where we used the optimality of 2d and of zhn We conclude that
2 { v (Er))
1 ( v. (2x))= R(Exy + An20)) + + A, ))) &
Taking the limit establishes 18 Since z" was an arbitrary element in a Āll, we obtain that the limit of A is
argmin (z, (v2 2)
subject to z 0 |A]]
where we used 16 Finally, observe that z" was an arbitrary limit point of zd and showed that all limit points are the same; this is enough to conclude the result.
D PROOF OF PROPOSITION 6
The proof of this statement relies strong convexity of Fn. Similar results have been proven in the context of entropic optimal transport (e.g. Genevay et al. (2019 ; Mena & Niles Weed (2019).
Our proof is similar the approach taken in Rigollet & Stromme 2022
17

---

Published as a conference paper at ICLR 2024
Let foo goo) minimize
Hoo
Note that Poo with
2
(DA fxolx) + (v0(11)
Poo = exp
minimizes
i0T (0
Let
Poo
n² yyoo.)
Poo (x)o()))) and Go = (goo(y))))
Note that by optimality of Aco, l | this be seen by comparing the objective (iT (ft) Aco and A. Moreover, due to the uniform bounds on &oo &cx from Lemma 12 Poo uniformly bounded away from 0 by exp(-C/ for some constant C that depends on â.
Let Pn minimise
n we know is of the form
2 An F Gn)
Pa = n²
exp
for vectors An, Fn, Gn- The 'certificates' are ooo and Zn pn where
Poo
Poo (poc and Φn 7 (Pn -Pn)
Note that 60 Px VA .t). The goal is to bound ||00 yo 1 4 n that nonde-
generacy of Zx would imply nondegeneracy of p Note that o Pn = s nn So, by the triangle inequality,
|00 p:0 p0
pn $ 1100 (poo @ () - Paloo + 1100 (fn - (f)-0
D 10 - Pn b + 1 1100 (p - Poollo + + 1100 (ńn f -1
The last two terms on the RHS can be controlled using Proposition 20 and are bounded by @(n -12) with probability at least 1 ((ep((- for f >> 0. For the first term on the RHS, letting Z Poo Pn.
| 10.0
Σ 1.
C(x, C(x,,, y,)
< ||Clloo
1
exp
( An Fy Gn)
exp (2 Foo @ Goo)
n f.j
Let _ ,( + Fn Goo vI + Foo @ Goo)lloo). According to the Lipschitz continuity of
the exponential
2 A0, A00 Foo @ Goo) i.j
exp
An + F Gn))
-exp E
exp
(19)
[.]
eexp(L/e) en
( An P G (0, Aoo + Poo + Gool)
(20)
12exp(L/e)
-n Aook + - _ -ood - I Good?
(21)
1.j
n
By strong convexity properties of In and hence Hn, it can be shown (see Prop 14 that 21 is upper bounded up to constant by
exxp(Lle) (Fo, Goo Abo) Kr Gn. An))
exp(LIe)
1
| 10 Gou Aboo)l?
4
18

---

Published as a conference paper at ICLR 2024
Inl Foo, Goo) |2 1
dā In Goo) 1²]
where &oo- = - (p B1) E A By Lemma 11 and Lemma 12 L 6( with
probability at least 1- (
in Finally,
Goo) + oz)
| Aocol 1
<2 80 al| 2 + Po - (Po /|}
1 1
1
R Σ j=1
n 01 In Foor Goo)
_
Poo yj)
11 1=1
n Fn Foo Good 1/
=
1
Σ
1
-
-
R n
Poo (xi.
19 and that these are bounded by Ø 1) with probability at least 1-0 (-((exp(--x)) and from Proposition 23 assuming 2 log(2s), (I 2,, S with probability at least 1 - 0 (e)p(-n)).
show Propositions
So, for some constant C> 0, A -0 - 1n S exp(C||²l) (e) with probability at least 9x-
1-exp(-r"), The second statement follows by combining our bound for 21 with the fact that
[nn ,,;,),,,,,,0111
In the following subsections, we complete the proof by establishing the required strong convexity properties of In and bound V fn 2 (Aze Fos Goo) using concentration inequalities. The proofs to some of these results are verbatim to the results of Rigollet & Stromme (2022) for deriving sampling complexity bounds in eOT, although they are included due to the difference in our setup.
D.1 STRONG CONVEXITY PROPERTIES OF Fn
In this section, we present some of the key properties of fn. Recall that since we assume that B are compactly supported, up to a rescaling the space, we assume without loss of generality that for all k, Cx(x,y) 1.
Lemma 11. Let A, minimize 1OT Assume that for some constant C1 > 0,
texp(C]/E)n 2 Amin [1, All))
then with probability at least 1 1 (exxp((-rr))
AA = ||||A||
Proof. Let A be minimizer to (0 By optimality of A, Аl L(A,A,) & A f)) Writing P,, _ ld and 0 [ A,Cg(x,.,y))
"
A A|l - (P,, A) (P, 0, H(P) & A| (A,,
(22)
2
for all P satisfying the marginal constraints P1 = 1 and PT1= 11.
Note that since (A. + K
-
(23)
=- f f,) Wa, (A) (a( KL @ ().
(24)
The first two terms can be shown to be (n -b Indeed.
R where Z A (iA, C(x),))) - EKÇĀ, ((x,))1I)
(0 - = Σ Z4
1=1
the sum of terms with mean zero. Moreover, |2,1 2 ||] By Lemma 18 ($7, át an) = with probability least -ex
/m
19

---

Published as a conference paper at ICLR 2024
The bound
We, (A) - W+(A) = =0(exp(-C/eem
with probability (exp((2)) is a due to the sample complexity of eOT Rigollet & Stromme 2002 obtain
Plugging this back into
22
A|l = A | On #)+ (n++))+
(H(P)- @ ()) + ((P,. 2n + A2)
(25)
remains to bound the terms Ti and Tz. Intuitively, if = pa e f has density p. then we can show that 7,, T2 = =( $) by choosing P Qn è (y)x,, However, Qn will only approximately satisfy the marginal constraints Qn = _ and Q = 1 (this approximation can be made precise using Proposition 16 So, we insert into the above inequality p= Qn with Qn being the projection of Qn onto the constraint set ( ( n. " n)).
By Altschuler et al. 2017 Lemma 7), the projection On satisfies
-1
Q 1
1| Qxlli = Qn 1---11
+2
(26)
By Proposition 16 with probability at least - 0 (exp(-xx)) 12)),
0n
1 |,
1 - n
ト _
p(x,y))
pix,,
=
n
11
n
Similarly, 1Q
n 11,,, On ż) with high probability.
Note that
n
k - Qn. pnAi| Σ Σ (nC Umli.
(27)
k=11,/=1
max IAAl Qn Qull
IAAl
(28)
Moreover,
( - Qn, An & |A| 0 Qn)
(29)
By Proposition 19 and 20 with probability at least 1 - -(xexp(-1)))
E| P, - d6 = 1) and E| - - = O 1)
So, we have T2 Cn ||l
We now consider the term Ti in 25 Since P is uniformly bounded from above by exp(C/e) and from below with constant exp(-C7e) for some C O, one can check from the projection procedure Altschuee et al. 2017 that Q, is also uniformly bounded from above by _ exp(C/e) and away from zero by n exp(-Cle) for some C So,
ITi & H((2)) H(Q.)) H(Q KL. -
@ po
(30)
e -.l. - B)
(31)
We can use the bound 26 to see that || - Qnll 2 t/v with probability at least 1 exp(-13) To bound T3 A H(Q,)) KL a @ 6), note that Moreover,
log(p (1 hatp://nnnkuonnno//ppoounhhhaaaa -ogg///t..//)))btlaaa.l..
E(Ta]
Σ
n²
1=1/=1
n(nn 1)
1)
Elog(pptz,,)))))(L
n²
20

---

Published as a conference paper at ICLR 2024
nº / log (ldcx 1 log(p(x, da
n/2
1 n (1pl.0 - 1) (o
Therefore, by the bounded differences lemmal with
n n
f( 2n) Σ
log(p(xx.yy))Dtx...)>
where (z). Then, letting 2j = = for all j # i, the bounded differences property is satisfied
with
f(2)..., f(z(
(, 1og
n²
i=1
-=1
4(
A[1pplls + +)
- c.
n
So, with probability least 1 2exp (-t)), 17al # In summary, with probability 1-0exp((-11000
2011 +1)
(A C2texp(Gi/e)n #) [[ & >||
(2(exp(C)/e1
vn
for some constants C₁, C2 > 0. So, choosing
C2texp(G//e) min (A/4,A All, /2).
we have
|| 2 A
with probability at least exp(--r0).
Lemma 12. Let (A, minimize Kn) Let = = Then,
F [-3 Callo and G E Calll
Moreover, exp(4[[CAllo0
> exp + (G+ > exp(-61lCalloe Note that
||C & |All
If(A.f.g) minimize
Let = AA Then, for x,y Koo
[ and -llellllllll Moreover, exp(4 ) exp (f (fe+g+ e
Proof. This proof is nearly identical to (Rigollet & Stromme 2022 Prop. 10): Let A, F.G minimize
[.
By the marginal constraints for P, na
(exp (ex given in p R
2 (o Fi Fl+
1= n
Σ
exp
(32)
2 (-1laall + +Fi) - exp (2G, = eep + +F:)
2
exp
-
21

---

Published as a conference paper at ICLR 2024
where use Jensen's inequality and the assumption that I, Gj 0 for the second. So, Pi rl for all E [n]. Using the other marginal constraint for Pn along with this bound on P. implies that
1= 1 I exp 11
o ()
exp (211cAlo++ + GD1)
+G)]
So, Gn,j > -211CAlloo
To prove the reverse bounds, by duality between 9
show that Σ; Fj is lower bounded: Note that since H(P) -1,
and
An
||CAlloo + (* P)
E = n 7
H(P)=
-
& E exp (F: GT+ [
2 ((
2
2n 1.
By assumption, Σ Gj 00 0 and exp (1F; ; Gj + (0, A) n²² So,
1 n Σー IIC
By repeating the argument in 32
we see that
a exp(/
+ Gf))
Σ Fi
> =cxp(2/((221Calle GJ)).
Ne
So, Gj CalCal and Fj -3
The proof for Koo is nearly identical and hence omitted.
Similarly to Rigollet & Stromme 2022 Lemma 11), we derive the following stroonngiggssongggggoo convexity bound for Fn:
Lemma 13. The functional In strongly convex with
f,(A,,.G)> IVA & (,,, G)),
exp(-L/E)
GTB²)
(33)
n 0n
E
n
n
for some -= v l.
Proof. To establish the strong inweetitiitiiuuuulllll let
h(t) A ( 0) + 1)
1g")
suffices find 8 such that for all (0,1),
((()) 6 1
||
(34)
Let Z, A ((1 ()F ( (G') (a1 (10)0)++00+0000 1DA). Note that
2
2
h²(t)
Σ
exp
(35)
E I.j
Since llcalloo v Vllcallo llo & L, by Lemmal
12
11 Sllcallo v 1lo So.
2 k1()) En
(xp((Lle))
Σ (F F; G G. + (0 A0) A)
1.j
By expanding out the brackets and using Σ, G; = and since C. are centred (C(C)), = and [C(C)), 0), (34 holds with 8 = exp(-L/e).
Based this strong convexity result, have the following bound.
22

---

Published as a conference paper at ICLR 2024
Proposition 14. Let (A Gn) minimise Fn. then for all ( of such that exp + G (, satisfy the marginal constraints of n - Fn. Gn)
L)p(
1
11
An)ll
+
IF
Full2
-| G2ll
e
122
n n
and
(r(A,, (An Gn) eexp(L/e)
+n 1100:33/144F0033++1//.3311111111
where L AI] 11)
Proof. By strong convexity of Jn. we can show that for any (A.f.g). (f' f'g') E y and any EEdlAll
F',
V (A.EG) (ALG))AAFFFFF.CCC G*) 4
ACC, A' A)
1) || - G 3
(36)
where = 2expl Lic) with I. = A(IA v [|])) The first statement follows by letting (A,FGG)
(An Fn, Gn) in the above inequality.
To prove the second statement, let
M CA + A'- A' A)
|x(AA - 1312 1 n
1 1 -
By minimising over (A', ("',F, note that
1 (r/1 [1 +I 26
M
So,
1
- M &
n ||
28
Finally, note that - = (.A.F. Kn(An, Gn) by optimality of An, Fn.
D.2 CONCENTRATION BOUNDS
Lemma 15 (McDiarmid's inequality). Let f f R satisfy the bounded differences property:
sup [x i fx)....,, 1, 1, i (,l,,.,,x))) xn) s. x/e3
Then, given Xi....,,; Xn random variables with X; E 2, for any t> 0,
Given random vectors and Y, denote Cov(X, Proposition 16. Let n have marginals and and let p marginals a, B. Assume that || pl & b. Then,
= =)
EE[X]> and Var(X) Cov(X, X). Let (t) ft where ft has
dx dapp
(b+1)
p(xn
n
n
and
1 n
Ø Σ p(x),)1)
(I+ b+ 1)2
Σ
n
j=1
n
n
with probability at least -
121(1)1))
23

---

Published as a conference paper at ICLR 2024
Remark 17. The bounds also translate to an & norm the marginal errors, since by Cauchy- Schwarz,
E p(x,
Σ
1
p(x
1
n1
Moreover, by Jensen's inequality E - 1 - y & /
EE 1 (1- - p(x)))))
Proof.
11
22
n Σ
p y7)
n²
ij.k=1
For each jE [n], we have the following cases for uj (11 (1-pl ()
i= == then " = - 1)2. There such term.
2. i= and k * j, then 2y = E(x,, (1- y1) p(z, p(z, O. There are such terms.
3. i + and k j, then U = E(e,y) (1- y)) p(z, y)) 0. There are n such terms.
4. i and # then uj = =Ex-a,y-p11- plx, y))² and there are 1) such terms.
5. i.j. 1 all distinct. Then, 2j 0 and there are (n - 1)(n - 2) such terms.
Therefore,
Σ -(00-pxx,)))
E 112 n²
-pxx Ix-ayyp(1-- p(x, y))
n
1j.=1
Using | p(x, y & b+ gives the first inequality.
Note also that letting V , - (11 p(x) p(x. yyD) 11 - " (1 (-x(xx,)) We will
apply Lemma 15 f(,....., 2,) = [V] Let = = - (1-p(x 1DD), where xi. = x for
i,j 2. Then, for vectors u of norm 1,
(u, V' - V)
U Σ
(p(x ,,,,
i=1
n1
Σ (p(x),>i) -1(x;x)()))))
Σ uj m (p pp(x),)))
j=2
2b
2b
Σ
Uj
Vn
So, by the reverse triangle inequality, |II v'll 2b Vn
It follows that
n 1/2 VII sn 1/2 1 1/2
EIVI 1/2 1+ -EIV 112 b+
with probability at least (x)))(= as required.
24

---

Published as a conference paper at ICLR 2024
D.2.1 BOUNDS FOR THE COST PARAMETRIZATION
In the following two propositions (Prop 20 and we assume that 4n is defined with Cg- (x, yp)), and discuss in the remark afterward how to account for the fact that our cost C₃ in are centered.
Note that the following classical result is direct consequence of Lemma] 15
Lemma 18. Suppose ZA.... Zm are independent mean zero random variables taking values in Hilbert space H. Suppose there some C' such that for all 7 C. Then, for all with probability at least 2exp(-1).
m
2
CC
m
k=1
m
Proposition 19. Let C maxx,y C|,,, y0|| Let (xi, (1) iid drawn from Then,
VBCt
| -00 x0"
n
with probability at least -exp(-x")
Proof. Direct consequence of Lemma 18
Proposition 20. Let 0. Let (xi, Vi be ii.d drawn from ât, which has marginals a.B. Let
P= (p(x. wheren has marginals a. B and dx
Taaø Then, El P (pa B) 1²
6n and
1+ Vn
|0 の- pa @ Bl|
with probability at least 1 2ex ||
Proof. Let h(x, y) & yCCx,,), then
$ P
h(x,yy)) Σ h(x;)1))
and
E(0" pa B] AA' (pa pa p)
n
follows that
E E|00; 0* pa pa
1
40 (pa pa B]] 2 Var(O;,))
n²
and
g Σ i.1.
Var(O), P)
yo)(h( yk)
n
Note that Cov(
whenever k, are distinct and there are n(n 1) 6n such terms, i.e. there are 6n² 12n + 6n nonzero terms in the [[C]] and are uniformly bounded.
2n 3) -6 12n sum. It follows that
To conclude, we apply Lemmal with f(z)..., -n) || P 4b pa bl| and z; (xi, v:). Let X ( 2n) " th pa & p. Then, Note that for an arbitrary vector u,
(u, X(x), X(Z),,2........ 2e...., 2n) X(x), 22. ................................ mnd)
Σ
Σ
us (h(xi, h(x/, 17))
[.]
Σ "k (h(x -h(x,) +
Σ
ux(h(x,,)) - n(x).,991)
k=
25

---

Published as a conference paper at ICLR 2024
< 88 ulll sup ||[xx, y) |p
x.
So, f(2,,22 2g....., 2n) f(2 22..... 2n) = 8n || and
1+t vn
|0 x pa
Bll -
with probability at least 1 (6x
Remark 21 (Adjusting for the centred cost parametrization). In the above propositions, compare " P= E Cy(x),y))P with o' x (x, y) y) for being a discretized version
of π. The delicate issue is that in particular,
Pn
Cg is centralized version of C₂ - (C,(x,y))1)) In
1
C C= C.
(CR).
Σ
(CWi.
ΣΣ (CD))
n
n Fli=1
Note that if P1 1 and p 1 then n
n
E(C = I(C.),)PI.) ni F
(Cx)..)
[.J
This last term on the RHS is negligible because @" (a @ (3 = = by assumption of C. being centred: Note that
t1 1
L Σ (CWR.
n 1xn)
nº
1
Applying the above proposition,
m+gog(n)
レ j=1k=1
(CD)k -
S
n
with probability least 1- ex) So, up constants, the above two propositions can be applied even for our centralized cost C.
D.2.2 INVERTIBILITY BOUND
Recall our assumption that M A -ux y(x, y) is invertible. In Proposition 23 we bound the deviation of D, from M in the spectral norm, and hence establish that it is invertible with high probability.
We will make use of the following matrix Bernstein inequality.
Theorem 22 (Matrix Bernstein). Tropp et al
2015 Let Zi...., Zn Rxd be independent symmerric & I. for all i [n]. Then, for all t = 0,
mean-zero random matrices such that 7
1 2olog((dd + Llog(2d)
E
Σ
2
where 02 = =||>|, Proposition 23. Assume that + +1 n. Let 1 00. 0. Then,
F|2
_ nス
1
Pn M 5
log(2s) f
n-1 1 vn
with probability at least ((exp(-12) 12)).
Proof. Recall that A EE where Ct is the centred version of the matrix č (E)(x,,)))ij.- That is,
(Cx),,) = (Cx)i.)
- 1 (CC)).m 1 [I(C))
[
(37)
n 7 1
n²
26

---

Published as a conference paper at ICLR 2024
For simplicity, we first do the proof for (L, =Ix AgCg(())))))))) and explain the end how to modify the proof to account for 4 using the centered C₃. First,
n² 0n - =
1
Σ
R1
Σ
C(x,
1, - C(x,, )) dp(r)
(38)
11 i=1
+ cl,, aply) c(x,)) - da(x)dp())
At
Σ
(39)
To bound the two terms in 339 let Z _ (Cx,)))((x, ( - c ()C,, y) da())
and observe that these are i.i.d. matrices with zero mean. By matrix Bernstein with the bounds || & and | 44,
8log(2s)
2log(2s)
log(2s)
E
Z
4
3n
assuming that log(2s) E n.
For the two terms in 38
1 C(x, C(x
C(x,,) ddps))
E
(40)
1
C(x., J.C(..)))
+E
C(x,, y)
Cixi, JICKK, ddpy))
2
(41)
2
11 +
Clx,, C(xx,)))((x.) - C(x, TT adpr)) jei
-
(42)
n2
For each = I, let Yj= yyxx, faari)y-CCCxxx) C(x, y) TT dp(y) and observe that condi- tional on xi, {Y1} jæz are iid matrices with zero mean. The matrix Bernstein inequality applied to R-T E jemnin Y, implies that
8log(2s)
2log(2s)
log(s)
E
Σ
& 41
n
jE(n]
Vn-1
3n-3
assuming that log(2s)
Finally, we apply Lemma 15 to (n))...,,2,) = , -
MM. - Let 0',u= Ex C,C,), with
x/ = x; and y/= y; for For each vector of unit norm,
1 (()) 2 - (00)) Du u) = n
1 n
yj) TT u1(² C(xk, up²
(43)
C(x;,))) T u12 - y/) T ul2
Σ ((),, T - up" < 4n (44)
+ n
So,
log(2s)
n 1 dD
M
<8
n1
with probability at least -exp(-2r2/16)
27

---

Published as a conference paper at ICLR 2024
To conclude this proof, now discuss how to modify the above proof in the case of the centered cost Ce, that is 4, Lx AxCx where Ce is as defined in (37 Note that in this case,
k
1 n
Σ
Σ
Ce()
Σ Ce(xi.y8))
(45)
na
1,f1
p..=1
n C2(xx,
n
n
Σ
Ce(Xp, vj)
Σ
Σ
Colx
Σ
Cg
n
n
1=1
(46)
We already know from the previous arguments that
Cglxx, M 0 llggsssmm 1,f=1
For the last term in 45
let A ji.j.p,4) q): ie [x] [ ( (m)
Note that A has n 1) ((((222))0 2n 3) terms, and A [m] m x [m) x [n] has (n²)) terms. Therefore, we can write
Σ
Clxp. Yay
Σ y
p.9=1
1,1=1
C(x,,))) T (i.j..q
Σ
+E
Σ (i.i.p.)cc-
Ca
n¹
n²
Σ
+0(n')
n²
(.J.P.Q)EA
where the second inequality is because there are (nz2) terms in A and we used the bound that |Cx, v) = Moreover,
Σ
C(xx.)y))Cxxx)))))
n=3 Σ Σ
E
Σ n- qE(L.).p)
1
C(xp. Coxx.1)))(x))
UIJ.p.qieA
n
JEUI
Note that conditional L.
Σ
C()x),))) 1= (lep..)d4(()) - =0
qE(1.).P)
by assumption that C(x, y)yp(y) O. So, apply Matrix Bernstein as before to show that
El
Σ
C(xx. 1S log(s)
n
(.J.P.qIEA
A similar argument can be applied to handle the two terms in 46) and so,
n 0 M S log(s)/n.
The high probability bound can now be derived using Lemma 15 as before.
E PROOFS FOR THE GAUSSIAN SETTING
Simplified problem To ease the computations, we will compute the Hessian of the following function (corresponding to the special case where Σ, = Id and Σ. = Id):
E
W(A) (Ap E) logdet(ld-
(47)
To retrieve the Hessian of the original function note that, since logdet(Z - Σー - Σ) logdet(Ep))
logdet(Id = I H
- B } a change of variable 2 x xΣ & shows that W(A) W A
and, hence,
v = (2) & A 2
(48)
28

---

Published as a conference paper at ICLR 2024
By the envelope theorem, (W(A) Σ, where Σ is the maximizer of 47 tation V² (A) to differentiating the optimality condition of Σ, i.e.
thus reducing the compu-
A= -E-1E(d
(49)
Recall also that such a has an explicit formula given in 5
Lemma 24.
Ww = ( = (ld (E--EE)--1 (AT A)T)
(50)
and Σ is in 49 and is the linear map defined by vee(X) ve
of Lemma 24 Define A) & E(d - Note that a maximizer of 47 satisfies G(,,A) = 0. Moreover, dgG is invertible at such (E, because this is the Hessian of logdet(Id= Σ 2), strictly concave function. By the implicit function theorem, there exists
function f: A Σ such that GAf(A) = and
ff = (0=G) -
It remains to compute GyG (f) A):
At an arbitrary point (E A) we have
G (Id- - Σ) ld (Id & 2)d = (ld - Σン @ ld (Id (ld-- 2) xx Idd (Id+
and at (f(A), A), i.e, with E(Id - A, we can further simplify
d=G= (ld Σ -ーΣ (Id = Σー
ATT
By the Woodbury matrix formula,
(Id E(Id -[Id=-" (Id-EE - =Id AΣT.
So,
- 1 (D I 1+ (AT $ AT,
thus concluding the proof.
We remark that from the connection between ((A) and (A)), i.e., 48 we obtain Lemma 7 as corollary.
E.1 LIMIT CASES
SVD representation of the covariance To derive the limiting expressions, we make an ob- servation on the singular value decomposition of Let the singular value decomposition of A be A= = TDVT, where D is the diagonal matrix with positive entries d,. Note that A = (Id+ (AT A
V(Id 2 Moreover, and A A commute,
22 D12 == == ADLEAA _ ATT = =U ld
- =nn VT =UŌVT,
(51)
4
where D is the diagonal matrix with diagonal entries
2
d, =
+
(52)
4d
2 d
29

---

Published as a conference paper at ICLR 2024
E.1.1 LINK WITH LASSO: LIMIT AS
Note that
d
E
1+
4d
d:
d +E
-,
-
(53)
2dj
E²
2d
E
E3
It follows that lim 1 0 and hence, eV WA) ld. So, the certificate converges to
E (2 Ealu.n) sign(A)). =
(54)
Proof of Proposition The iOT problem approaches a Lasso problem as E Recall that in the Gaussian setting, the iOT problem is of the form
argmin & (A) ||ll (A, Σe,A 2. 2
logdet (xpleetIp-piieeee I Ee,A)
(55)
where Σe, satisfies ΣB
-T E,
- Σ = E2 -E,1 1 So, we can write
argmin a _ 1 + +A, E,, - - &
2c,1A-1
(56)
A
Let
x イ and =E E
From 52
X has SVD decomposition W U diag(d:) then
Σ=,A= where D = diag(d;)
and
d d E
+Gle 5
So,
1ogdet(<2g Ee,1A 1 ) logdet E SE,1 1
logdet (U diag(1
dil +0 -
= logdet (diag 1- -0(
= = 1XXI (1)
O0E
Also,
E(E,,A -2 L,, A) E X) (X Xo. x)+
So, assuming that 1 11=000 Ag/E,
(F) = 11 |All
(X - X X) || + +2
(57)
= A
(2) の 0 (A- A)
- 12) 2 A +
(58)
The final statement on the convergence of minimizers follows by Gamma--onvergence.
E.1.2 LINK WITH GRAPHICAL LASSO
In the special case where the covariances are the identity (E, = Id and Σ3 Id) and A symmetric positive definite have that Σ is also positive definite and Galichon's formula (5) holds (since A invertible) and hence the Hessian reduces to
- 42 (wW))
(A A) & 3 +T.
(59)
Moreover, if A admits eigenvalue decomposition = UDUT, then Σ admits an eigenvalue decomposition = T with entries of D given by (22) Note that it follows that lim,-0 - = Id and, hence, lim (Σ & Σ Id T, singular matrix with the kernel being the set of
30

---

Published as a conference paper at ICLR 2024
asymmetric matrices. So, the limit does not necessarily exist as 0. However, in the special case where is symmetric positive definite, one can show that the certificates remain well defined as - O: Let S- be the set of matrices such that y is symmetric and S_ be the set of matrices y such that w is anti-symmetric. Then, (V² ((A) (S=)) c S+ and (V2)(A) *(S_) c S-. Moreover,
Is. (A@ () A)
3o +Id).
()w))))) Is- (A - -
Since the symmetry of implies the symmetry of sign(A), replace the Hessian given in 59 by (A 8 A) S 1 +ld and, hence, since lim -= Id, the limit as E 0 is
lim & = (A- @ A Jennaa a..)) sign(A)). 00
(60)
This coincides precisely with the certificate of the graphical lasso:
argmin θ) -oogdet())) + 11101111
6z:
Proof of Proposition The iOT problem with identity covariances restricted to the set of positive semi-definite matrices has the form
arg mig Fr, (A), where Fen( & Alllh + (A, Σ logdet(ld Ze,
(61)
A=0
where 1-Σ7 - EL From the singular value decomposition of L....... 51 that if A symmetric positive definite, then so is Σg,s- Plugging the optimality condition, i.e., T ZE, A into we obtain
E argmin | + (A, Σe,A - E) logdetleA EE,A)
AzD
(ogd (Ee,A - Σ)1E) + E logdet((e,.))
argminnn |AIl
AzD
2
Let A = Edo for some Ag 0. Then, removing the constant (logde((eld) term and factoring out the problem is equivalent to
arg ||l logdet(A) (A, '(E,A Σ)) logdet(EE,.) A=0
Assume that f -- From the expression for the singular values Σ.,1 in note that
Σg Id -a--dd--+---e +0-2+)(e3)
So, lime-o((E,,' - Eg,i)l== -2, A-1). - (A-1-AAp
The objective converges pointwise to logdet(A) A) - A-A-
Ao |1Al
and the statement is then a direct consequence of Gamma Caamma-convegeene convergence.
Remark 25. Note that from
52 we have I= =Ild & (EE2). So, in this case, the covariance of
is
ld
Id - A
ld Au-5AA11 - (E
ld
For (X, Y) π, The Schur complement of this is the covariance of X conditional on Y, which EA-i + +E)).. So, up the e, one can see -1 as the "precision" matrix the covariance of X conditional on Y.
31

---

Published as a conference paper at ICLR 2024
F LARGE SCALE -IT SOLVER
Recall that the iOT optimization problem, recast over the dual potentials for empirical measures, reads
1 +r (0, 2n 1.j
& Σ exp
2 (Fi + Gy+ (0 +
inf
A.E.
To obtain a betterr-oondiiooned optimization problem, in line with/Cuturi & Peyré 2016 we instead consider the semi-dual problem, which is derived by leveraging the closed-form expression for the optimal G, given F.
inf A.F
E P: + (RnA) + - Elog log
n Σ
- (F: exp
((FF)(( 1
Following Poon & Peyré 2021), which proposes state-of-the-arr Lasso solver, the last step is to use the following Hadamard product over- pareeettatttiiniionnnn of the e' norm
111/
II 2
All
min UeV
where the Hadamard iroduct o V & (U, Vi),, to obtain the final optimization problem
log Σ exp (Fi (U inf 1 S F: (0, V & _ A,U.V
2
(U, VMn.)) + uuns || V
This is a smooth optimization problem, for which we employ a quasi-Newtor solver (L-BFGS). Although is non-convex, demonstrated in Poon & Peyré (021), the non-convexity is benign, ensuring the solver always converges to a globba minimizer, this, one can reconstruct the cost parameter, A & U ④ V*.
U V*), of the functional. From
32


Reviews:
[{"review_id": "ssoAub2Ktg", "reviewer": "Reviewer_gfhu", "summary": "The paper is about inverse optimal transport (iOT), that is the task of inferring a ground cost from the sampling of an (entropy-regularized) optimal transport plan.\nFor linear parametrizations of the cose, iOT is a convex inverse problem (Dupuy et al 2019), that the authors propose to regularize with the $\\ell_1$ penalty.\n\nThe paper first derives an irrepresentability condition (IC) for the iOT problem,\nbased on non degeneracy of a \"precertificate\" (Def 1).\nFirst, the authors show that the solution of the regularized \"infinite samples\" problem shares the same sign as the true parameters under IC, for low enough regularization value (Thm 3).\nThen the more interesting \"Sparsistency\" results are then derived: in the finite sample case, the correct support is still recovered (under IC) for small regularization values and sufficiently many samples (Theorem 7).\n\nAdditional details in the case of Gaussian distributions are provided, with interpretation as Lasso and Graphical Lasso respectively for vanishing or exploding entropic regularization strength.\nExperiments on limited size graphs (80 nodes) conclude this theoretical paper.", "questions": "Questions:\n- why does Proposition 8 imply that $z_\\infty$ is nondegenerate? Why, given its definition 1 line above, is it of the form written in (PrC)?\n- In prop 9 shouldn't the limit of $\\epsilon_n$ be $\\infty$?\n\n\nMinor:\n- In theorem 3, $\\hat A$ is such that $\\hat \\pi = \\text{Sink}(\\hat A, \\epsilon$; this has been introduced earlier but given the large number of variables defined in the paper, it may not hurt to recall it here; in addition, before, Sink was applied to $c$, so it should be $\\text{Sink}(c_{\\hat A}, \\epsilon$. In prop 10, prop 11 the order of the arguments is reversed.\n- can the authors explain the name \"precertificate\" (what's \"pre\" about it? Is there a difference with what Dunner and coauthors call \"Dual certificate\" in https://arxiv.org/abs/1602.05205)?\n- The authors refer to \"condition PrC\", but that is just an equality. Condition PrC would be that the dual norm of the precertificate outside the support is $< \\lambda$\n- what's \"the convex dual of W (A)\"? convex conjugate of a function (as in Proposition 4)/convex dual of a convex optimization problem?\n- is the notation $\\langle c, A \\rangle$ for a non scalar result?\n- In theorem 7 I spent quite some time looking for what the number $m$ was in other parts of the paper (because of the analogy $m/n$ in compressed sensing) before realizing that it was a free paramert; maybe replacing it by $\\ln(1/\\delta)$ is more common (that is only a suggestion).\n- below 9, in $A_n$ definition, both $\\epsilon$'s should be $\\epsilon_n$? Same in Prop 11\n- it follows (see e.g. Hiriart-Urruty et al. (1993)) : can you point to a specific result in the book? Also in the bibliography, the authors names appear twice for this book.\n\nTypos:\n- the paper does not seem to use the unmodified iclr template, the font differs from that of other papers\n- in Problem iOT L1 hat Pi, the first argument of L should be $A$ not $c_A$ (see def of L 2 equations above)\n-  The iOT problem of iOT\n- in a series of paper*s*\n- Thm 3: the solution ... satisf*ies*\n-  minimial norm\n- sufficinetly$\n- result above implies that for provided that the number\n- as exposed in Section ??\n- Cauchy Schwartz\n- see Proposition PrC (it's an equation not a proposition, and I don't see why pRC shows that W is $C^2$)", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 1, "contribution": 3, "strengths": "Inverse optimal transport has recently attracted attention in the community due to its potential impact in ML.\nProposing better alternatives to solve this nonlinear inverse problem is thus of interest.\nThe paper provides results on both \"full distribution\" and finite sample problems.\nI did not spot any mathematical error, but could not check all of the paper.", "weaknesses": "The authors could really afford to improve the pedagogy of the paper, which is quite heavy in terms of notation. Theory is quite involved, require background references to other works such as Carlier or Galichon.\nThe experiment description is a big block which, in my opinion, does not bring as many insights as it could.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "bdIolQC1ef", "reviewer": "Reviewer_akqu", "summary": "The paper studies the regularized inverse entropic optimal transport problem. Inverse optimal transport (iOT) is the problem of recovering the ground cost given samples from the (potentially, noisy) joint distribution. Recent works have proposed solvers to solve the primal and dual formulations of the iOT and regularized iOT problem. This paper studies the recovery guarantees for the L1-regularized iOT problem. They further explore a special case where the densities are Gaussian, and show that in some cases, iOT results in the graphical LASSO problem.", "questions": "- How practical is inverse OT? I understand that the current work considers linear cost functions. Can one, albeit without strong theoretical guarantees, learn a general (perhaps, regular) cost function from the samples drawn from the coupling?", "limitations": "", "rating": 8, "confidence": 2, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "- Inverse optimal transport is an interesting problem and an interesting take on the metric learning problem, this work takes a significant step forward in establishing a theoretical grounding for the regularized iOT problem.\n- I really enjoyed reading the paper, the writing is very clear, the background, method, and the results are well presented.\n- Showing that graphical LASSO as a special case of inverse OT is very interesting and makes sense.", "weaknesses": "- Nothing that I can think of.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "vopjeVfbeA", "reviewer": "Reviewer_2C9j", "summary": "This paper deals with the problem of inverse optimal transport on compact (but not necessarily finite) state space to find the cost function from a given empirical probability coupling. It establishes a connection to lasso and produces the sample complexity of solving the corresponding primal-dual method.", "questions": "Can the authors please elaborate more on the penalty term $\\lambda$? Is $\\lambda$ given for a given problem? Is it shrinking with $n$? \n\nAlso, do the authors implicitly assume that we can sample from the coupled empirical density? Or do we just have access to some samples?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Finding the cost function from the empirical process induced by the optimal transport algorithm is an important open question. In that respect, this paper addresses an important open question. The connection of the non-degenerate precertificate assumption with irrepresentability assumption in lasso is an interesting insight. Finally, I found the connections to the SVD of the covariance matrix illuminating. The gaussian example also serves to demonstrate the theory through the lenses of an example.", "weaknesses": "Although the sample complexity bound is good, there was no discussion about the tightness of the said bound. I wonder if the $1/\\sqrt{n}$ is also the best lower bound. Observing that the empirical density acts like a ``plug-in\" for the unknown true coupled density, can the statistical guarantees of the plug-in be translated to the guarantees for the cost function?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "QOuJa13iLF", "reviewer": "Reviewer_rfzP", "summary": "This article addresses the problem of inverse optimal transport, which involves estimating the transport cost from an optimal (noisy) transport plan. The authors focus on the 'penalized' $\\ell_1$ cost formulation, parametrized by a matrix of parameters (the cost is parametrized as linear combination of individual costs, such as in the Mahalanobis setting). The article provides two types of guarantees:\n\nThe first type of guarantee, referred to as the \"finite sample case,\" is an estimation guarantee obtained from dual certificates when the observed transport plan is an entropic regularized plan between two empirical measures.\n\nThe second type, in the \"Gaussian case,\" deals with the scenario where the measures involved are Gaussian distributions (a well-known case in transport that admits a closed-form solution). In this case, the article demonstrates that the cost estimation can be solved as a graphical lasso problem, especially when the entropic regularization is small.", "questions": "I'm curious to know whether the Gaussian results really require entropic regularization. Without regularization we also have a closed form: can't we deduce good invOT estimation guarantees in the non-regularized case?", "limitations": "", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "- I find the introduction and contextualization of the article to be very well done. The related work appears comprehensive, and the problem is well-situated.\n\n- The theoretical results in this article are interesting. For the discrete case, the authors demonstrate that, with a small $\\ell_1$ regularization and a sufficiently large number of samples, the ill-posed problem of invOT can indeed recover the costs.\n\n- The connections between graphical lasso and invOT are also intriguing and interesting.", "weaknesses": "- I find that the article doesn't put in enough effort to explain the practical implications of the various theoretical results, which remain somewhat abstract. It's quite challenging to understand how these results can be practically applied.\n\nFirstly, the concept of a pre-certificate condition is rather very abstract (unlike the case of Lasso that we can somehow interpret). It would be interesting to provide the reader with a bit more guidance to offer a small interpretation of this quantity.\n\nAnother example is Proposition 8. It's difficult for me to extract something from this result, aside from the qualitative interpretation that \"for small regularization, under somewhat abstract assumptions of certificates, and with a sufficiently large number of samples, solving the invOT problem yields a good solution.\" To make these results more applicable in practice, it would be helpful to either provide experiments or guarantees that demonstrate how to understand these pre-certificate assumptions or ensure that the level of regularization is appropriate. I realize that these may be challenging questions, but I find that not much intuition is given, and the practical use of these results does not appear straightforward.\n\n- My main criticism concerns the experiments section. I find the experiments too limited and somewhat confusing, making it difficult to obtain meaningful information.\n\nFirst, only the Gaussian with an identity covariance is considered. The idea is to study the impact of entropic regularization on invOT estimation. The presentation is somewhat unclear, and I struggle to understand from Figure 1 how it quantifies the influence of $\\epsilon$ on the estimation. There is no legend for the y-axis, and while it's understood to be the certificate value between two nodes, it's unclear how it serves as a good performance measure for the invOT problem. \n\nFigure 2 is equally unenlightening. It aims to determine the influence of the number of samples on the estimation in the very simple case of a circular graph and Gaussian measures. However, what is the x-axis on these three figures? Is it the number of iterations? The geodesic distance? If it's the geodesic distance, it's not clear because the y-axis is a global performance measure over C, while the x-axis appears to be a local measure, so it's unclear how one evolves with the other.\n\nThe conclusion appears to be that the estimation is \"good\" when the entropic regularization is sufficiently large, and the number of samples is also large. More importantly, there seems to be a significant gap: when $\\epsilon \\leq 10$, the estimation is consistently poor. Is this a result expected by the theory regarding $\\epsilon$?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (poster)