# Sparsistency for Inverse Optimal Transport

Francisco Andrade
CNRS and ENS, PSL Université
<EMAIL_ADDRESS>Gabriel Peyré
CNRS and ENS, PSL Université
<EMAIL_ADDRESS>Clarice Poon
University of Warwick
<EMAIL_ADDRESS>
###### 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, we work out in
detail the Gaussian case. We show that as the entropic penalty varies, the iOT
problem interpolates between a graphical Lasso and a classical Lasso, thereby
establishing a connection between iOT and graph estimation, an important
problem in ML.

##  1 Introduction

Optimal transport has emerged as a key theoretical and numerical ingredient in
machine learning for performing learning over probability distributions. It
enables the comparison of probability distributions in a “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 of this OT approach to
ML is inherently tied to the 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 a priori highly
non-convex nature when framed as an optimization problem, thereby inheriting
complications in its mathematical analysis. As we illustrate in 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 a restrictive setup, but it arises in practice (refer to 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 to stably learn a
ground cost from partial matching observations.

###  1.1 Previous Works

##### Entropic Optimal Transport.

OT has been instrumental in defining and studying various procedures at the
core of many ML pipelines, such as bag-of-features matching [Rubner et al.,
2000], distances in NLP [Kusner et al., 2015], generative modeling [Arjovsky
et al., 2017], flow evolution for sampling [De Bortoli et al., 2021], and even
single-cell trajectory inference [Schiebinger et al., 2019]. We refer to the
monographs [Santambrogio, 2015] for detailed accounts on the theory of OT, and
[Peyré et al., 2019] for its computational aspects. Of primary importance to
our 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]. More importantly, it
defines a 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éonard, 2012]. The role of noise in matching
(with applications in economics) and its relation to entropic OT were advanced
in a series of papers by Galichon and collaborators [Galichon and Salanié,
2010, Dupuy and Galichon, 2014, Galichon and 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., 2013] 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 et al.,
2006, Davis and Dhillon, 2008].

##### OT Ground Metric Learning.

The problem of estimating the ground cost driving OT in a supervised manner
was first addressed by [Cuturi and Avis, 2014]. 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 and Guibas, 2012, Xu
et al., 2018, Kerdoncuff et al., 2021] or geodesic distances [Heitz et al.,
2021], to devise more efficient learning schemes. The study [Zen et al., 2014]
conducts ground metric learning and matrix factorization simultaneously,
finding applications in NLP [Huang et al., 2016]. It is noteworthy that ground
metric learning can also be linked to generative models through adversarial
training [Genevay et al., 2018] and to robust learning [Paty and 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 and Galichon, 2014] over a discrete space (also see
[Galichon and Salanié, 2022] for a more detailed analysis), making the
fundamental remark that the maximum likelihood estimator amounts to solving a
convex problem. The mathematical properties of the iOT problem for discrete
space (i.e., 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 as we focus on 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 a ill-posed problem, and [Dupuy et al., 2019, Carlier et
al., 2023] propose regularizing the maximum likelihood estimator with either a
low-rank (using a nuclear norm penalty) or a sparse prior (using an $\ell^{1}$
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 to [Ma et al., 2020] for an efficient
solver without inner loop calls to Sinkhorn’s algorithm. This approach was
further refined in [Stuart and Wolfram, 2020], deriving it from a Bayesian
interpretation, enabling the use of MCMC methods to sample the posterior
instead of optimizing a pointwise estimate (as we consider here). They also
propose parameterizing cost functions as geodesic distances on graphs (while
we consider only linear models to maintain convexity). An important
application of iOT to ML, explored by [Li et al., 2019], is to perform link
prediction by solving new OT problems once the cost has been estimated from
the observed couplings. Another category of ML application of iOT is
representation learning (learning embeddings of data into, e.g., an Euclidean
space) from pairwise interactions, as demonstrated by [Shi et al., 2023],
which recasts contrastive learning as a specific instance of iOT.

##### Inverse problems and model selection.

The iOT problem is formally a bilevel problem, as the observation model
necessitates solving an OT problem as an inner-level program [Colson et al.,
2005] – we refer to [Eisenberger et al., 2022] for a recent numerical
treatment of bilevel programming with entropic OT. The iOT can thus be
conceptualized as an “inverse optimization” problem [Zhang and Liu, 1996,
Ahuja and 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 is 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
a symmetric matrix for Mahalanobis norms). There exists a 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, 1996] (least square regularized by
$\ell^{1}$), which culminated in the theory of compressed sensing [Candès et
al., 2006]. These theoretical results are predicated on a so-called
“irrepresentability condition” [Zhao and Yu, 2006], which ensures the
stability of the support. While our analysis is 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 a vector) and the estimation process necessitates solving a
linear program of potentially infinite dimension (in the limit of a 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 of
regularized iOT estimation, focusing on the special case of sparse $\ell^{1}$
regularization (the $\ell^{1}$-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 6, which leverages this
abstract irrepresentability condition to ensure sparsistency of the
$\ell^{1}$-iOT method. This relates to the robust estimation of the cost and
its support in some linear model, assuming the number $n$ of samples is large
enough. Specifically, we demonstrate a sample complexity of $n^{-1/2}$. Our
subsequent sets of contributions are centered on the case of matching between
samples of Gaussian distributions. Herein, we illustrate in Lemma 8 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 9 and Proposition 10 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 5.2 further provide intuition about how
$\varepsilon$ and the geometry of the graph associated with a sparse cost
impact sparsistency. As a minor numerical contribution, we present in Appendix
D a large-scale $\ell^{1}$-iOT solver, implemented in JAX and distributed as
open-source software.

##  2 Inverse optimal transport

##### The forward problem

Given probability distributions $\alpha\in\mathcal{P}(\mathcal{X})$,
$\beta\in\mathcal{P}(\mathcal{Y})$ and cost function
$c:\mathcal{X}\times\mathcal{Y}\to\mathbb{R}$, the entropic optimal transport
problem seeks to compute a coupling density

$\mathrm{Sink}(c,\varepsilon)\triangleq\operatorname*{arg\,max}_{\pi\in\mathcal{U}(\alpha,\beta)}\langle
c,\,\pi\rangle-\frac{\varepsilon}{2}\mathrm{KL}(\pi|\alpha\otimes\beta)\quad\text{where}\quad\langle
c,\,\pi\rangle\triangleq\int c(x,y)d\pi(x,y),$ (1)

where $\mathrm{KL}(\pi|\xi)\triangleq\int\log(d\pi/d\xi)d\pi-\int d\pi$ and
$\mathcal{U}(\alpha,\beta)$ is the space of all probability measures
$\pi\in\mathcal{P}(\mathcal{X}\times\mathcal{Y})$ with marginals
$\alpha,\beta$.

##### The inverse problem

The inverse optimal transport problem seeks to recover the cost function $c$
given an approximation $\hat{\pi}_{n}$ of the probability coupling
$\hat{\pi}\triangleq\mathrm{Sink}(c,\varepsilon)$. A typical setting is an
empirical probability coupling
$\hat{\pi}_{n}=\frac{1}{n}\sum_{i=1}^{n}\delta_{(x_{i},y_{i})}$ where
$(x_{i},y_{i})_{i=1}^{n}\overset{iid}{\sim}\hat{\pi}$. See Section 2.2.

##### The Loss function

The iOT problem has been proposed and studied in a series of paper, see
Section 1.2. The approach is typically to consider some linear
parameterization of the cost $c_{A}(x,y)$ by some parameter
$A\in\mathbb{R}^{s}$. The key observation of [Dupuy et al., 2019] is that the
negative log-likelihood of $\hat{\pi}$ at parameter value $A$ is given by

$\mathcal{L}(A,\hat{\pi})\triangleq-\langle
c_{A},\,\hat{\pi}\rangle+W(A)\quad\text{where}\quad
W(A)\triangleq\sup_{\pi\in\mathcal{U}(\alpha,\beta)}\langle
c_{A},\,\pi\rangle-\frac{\varepsilon}{2}\mathrm{KL}(\pi|\alpha\otimes\beta).$

So, computing the maximum likelihood estimator $A$ for the cost corresponds to
minimizing the _convex_ ‘loss function’ $\mathcal{L}(A,\hat{\pi})$. We write
the parameterization as

$\Phi:A\in\mathbb{R}^{s}\mapsto c_{A}=\langle
A,\,\mathbf{C}\rangle=\sum_{j}A_{j}\mathbf{C}_{j},\quad\text{where}\quad\mathbf{C}_{j}\in\mathcal{C}(\mathcal{X}\times\mathcal{Y}).$

A relevant example is quadratic loss functions, so that for
$\mathcal{X}\subset\mathbb{R}^{d_{1}}$,
$\mathcal{Y}\subset\mathbb{R}^{d_{2}}$, given $A\in\mathbb{R}^{d_{1}\times
d_{2}}$, $c_{A}(x,y)=x^{\top}Ay$. In this case, $s=d_{1}d_{2}$ and for
$k=(i,j)\in[d_{1}]\times[d_{2}]$, $\mathbf{C}_{k}(x,y)=x_{i}y_{j}$.

#####  $\ell_{1}$-iOT

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 $\ell_{1}$-regularization [Carlier et al., 2023] to
recover sparse parametrizations:

$\operatorname*{arg\,min}_{A}\mathcal{F}(A),\quad\text{where}\quad\mathcal{F}(A)\triangleq\lambda\left\|A\right\|_{1}+\mathcal{L}(c_{A},\hat{\pi}).$
($\mathrm{iOT-}\ell_{1}(\hat{\pi})$)

##### Kantorovich formulation

By considering the convex dual of $W(A)$, one can show that
($\mathrm{iOT-}\ell_{1}(\hat{\pi})$) has the following equivalent Kantorovich
formulation [Carlier et al., 2023]:

$\operatorname*{arg\,min}_{A,f,g}\mathcal{K}(A,f,g),\quad\text{where}\quad\mathcal{K}(A,f,g)\triangleq\mathcal{J}(A,f,g)+\lambda\left\|A\right\|_{1},\quad\text{and}$
($\mathcal{K}_{\infty}$)
$\mathcal{J}(A,f,g)\triangleq-\int\left(f(x)+g(y)+\Phi
A(x,y)\right)d\hat{\pi}(x,y)+\frac{\varepsilon}{2}\int\exp\left(\frac{2(f(x)+g(y)+\Phi
A(x,y))}{\varepsilon}\right)d\alpha(x)d\beta(y).$

Based on this formulation, various algorithms have been proposed, including
alternating minimization with proximal updates [Carlier et al., 2023]. Section
D details a new large scale solver that we use for the numerical simulations.

###  2.1 Invariances and assumptions

######  Assumption 1.

We first assume that $\mathcal{X}$ and $\mathcal{Y}$ are compact.

Note that $\mathcal{J}$ has the translation invariance property that for any
constant function $u$, $\mathcal{J}(f+u,g-u,A)=\mathcal{J}(f,g,A)$, so to
remove this invariance, throughout, we restrict the optimization of
($\mathcal{K}_{\infty}$) to the set

$\mathcal{S}\triangleq\left\\{(A,f,g)\in\mathbb{R}^{s}\times
L^{2}(\alpha)\times L^{2}(\beta)\;;\;\int g(y)d\beta(y)=0\right\\}.$ (2)

Next, we make some assumptions on the cost to remove invariances in the iOT
problem.

######  Assumption 2 (Assumption on the cost).

  * (i)

$\mathbb{E}_{(x,y)\sim\alpha\otimes\beta}[\mathbf{C}(x,y)\mathbf{C}(x,y)^{\top}]\succeq\mathrm{Id}$
is invertible.

  * (ii)

$\left\|\mathbf{C}(x,y)\right\|\leqslant 1$ for $\alpha$ almost every $x$ and
$\beta$-almost every $y$.

  * (iii)

for all $k$, $\int\mathbf{C}_{k}(x,y)\mathrm{d}\alpha(x)=0$ for $\beta$-a.e.
$y$ and $\int\mathbf{C}_{k}(x,y)\mathrm{d}\beta(y)=0$ for $\alpha$-a.e. $x$.

Under these assumptions, it can be shown that iOT has a unique solution (see
remark after Proposition 2). Assumption 2 (i) is to ensure that $c_{A}=\Phi A$
is uniquely determined by $A$. Assumption 2 (ii) is without loss of
generality, since we assume that $\alpha,\beta$ 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 et al., 2020]:

$\mathcal{J}(A,f,g)=\mathcal{J}(A^{\prime},f^{\prime},g^{\prime})\iff
c_{A}+(f\oplus g)=c_{A^{\prime}}+(f^{\prime}\oplus g^{\prime}).$

As observed in [Carlier et al., 2023], any cost can be adjusted to fit this
assumption: one can define

$\tilde{\mathbf{C}}_{k}(x,y)=\mathbf{C}_{k}(x,y)-u_{k}(x)-v_{k}(y)$

where $u_{k}(x)=\int\mathbf{C}_{k}(x,y)d\beta(y)\quad\text{and}\quad
v_{k}(y)=\int\mathbf{C}_{k}(x,y)d\alpha(x)-\int\mathbf{C}_{k}(x,y)d\alpha(x)d\beta(y).$
Letting $\tilde{\Phi}A=\sum_{k}A_{k}\tilde{\mathbf{C}}_{k}$, we have
$\left(\left(f-\sum_{k}A_{k}u_{k}\right)\oplus\left(g-\sum_{k}A_{k}v_{k}\right)\right)+\tilde{\Phi}A=(f\oplus
g)+\Phi A.$ So optimization with the parametrization $\Phi$ is equivalent to
optimization with $\tilde{\Phi}$.

NB: For the quadratic cost $c_{k}(x,y)=x_{i}y_{j}$ for $k=(i,j)$, condition
(iii) corresponds to _recentering_ the data points, and taking $x\mapsto
x-\int xd\alpha(x)$ and $y\mapsto y-\int yd\beta(y)$. Condition (ii) holds if
$\left\|x\right\|\vee\left\|y\right\|\leqslant 1$ for $\alpha$-a.e. $x$ and
$\beta$-a.e. $y$. Condition (i) corresponds to invertibility of
$\mathbb{E}_{\alpha}[xx^{\top}]\otimes\mathbb{E}_{\beta}[yy^{\top}]$.

###  2.2 The finite sample problem

In practice, we are faced with the following finite sample problem: Given
$\hat{\pi}=\mathrm{Sink}(c_{\widehat{A}},\varepsilon)$, recover
${\widehat{A}}$ from the empirical measure
$\hat{\pi}_{n}=\frac{1}{n}\sum_{i=1}^{n}\delta_{x_{i},y_{i}}$ with
$(x_{i},y_{i})\overset{iid}{\sim}\hat{\pi}$. We will consider solutions to
$\mathrm{iOT-}\ell_{1}(\hat{\pi}_{n})$ and derived conditions under which the
support of ${\widehat{A}}$ is identifiable. The problem
$\mathrm{iOT-}\ell_{1}(\hat{\pi}_{n})$ can be formulated entirely in finite
dimensions as follows. Let

$H(P|Q)\triangleq\sum_{i,j}P_{i,j}(\log(P_{i,j}/Q_{i,j})-1).$

Note that $\hat{\pi}_{n}$ has marginals
$\hat{a}_{n}=\frac{1}{n}\sum_{i=1}^{n}\delta_{x_{i}}$ and
$\hat{b}_{n}=\frac{1}{n}\sum_{i=1}^{n}\delta_{y_{i}}$. We can interpret
$\hat{\pi}_{n}$ as the matrix $\hat{P}_{n}=\frac{1}{n}\mathrm{Id}_{n\times n}$
and the “noisy” primal problem can be equivalently written as

$\min_{A\in\mathbb{R}^{s}}\sup_{P\in\mathbb{R}^{n\times
n}_{+}}\lambda\left\|A\right\|_{1}+\langle\Phi_{n}A,\,P-\hat{P}_{n}\rangle-\frac{\varepsilon}{2}H(P)\quad
s.t.\quad P\mathds{1}=\frac{1}{n}\mathds{1}\quad\text{and}\quad
P^{\top}\mathds{1}=\frac{1}{n}\mathds{1},$ ($\mathcal{P}_{n}$)

where we write $H(P)\triangleq
H(P|\frac{1}{n^{2}}\mathds{1}\otimes\mathds{1})$ and
$\Phi_{n}A=\sum_{k}A_{k}C_{k},\quad\text{where}\quad
C_{k}=(\mathbf{C}_{k}(x_{i},y_{j}))_{i,j\in[n]}\in\mathbb{R}^{n\times n}.$
Note that the finite dimensional problem has the same invariances as
($\mathrm{iOT-}\ell_{1}(\hat{\pi})$), so, we will take $C_{k}$ to be centred
so that for all $i$, $\sum_{i}(C_{k})_{i,j}=0$ and for all $j$,
$\sum_{j}(C_{k})_{i,j}=0$. The finite sample Kantorovich formulation is

$\inf_{A,F,G\in\mathcal{S}_{n}}\mathcal{K}_{n}(A,F,G)\quad\text{where}\quad\mathcal{K}_{n}(A,F,G)\triangleq\mathcal{J}_{n}(A,F,G)+\lambda\left\|A\right\|_{1}\quad\text{and}$
($\mathcal{K}_{n}$)
$\mathcal{J}_{n}(A,F,G)\triangleq-\sum_{i,j}\left(F_{i}\oplus
G_{j}+(\Phi_{n}A)_{i,j}\right){\left(\hat{P}_{n}\right)_{i,j}}+\frac{\varepsilon}{2n^{2}}\sum_{i,j}\exp\left(\frac{2}{\varepsilon}(F_{i}+G_{j}+(\Phi_{n}A)_{i,j}\right),$

and we restrict the optimization of ($\mathcal{K}_{n}$) over
$\mathcal{S}_{n}\triangleq\left\\{(A,F,G)\in\mathbb{R}^{s}\times\mathbb{R}^{n}\times\mathbb{R}^{n}\;;\;\sum_{j}G_{j}=0\right\\}.$

##  3 The certificate for sparsistency

###  3.1 Certificate and non-degeneracy

In this section we present a sufficient condition for support recovery, that
we term _non-degeneracy of the certificate_. Under non-degeneracy of the
certificate we obtain support recovery as stated in Theorem 3, a known result
whose proof can be found _e.g._ in [Lee et al., 2015]. This condition can be
seen as a generalization of the celebrated _Lasso’s Irrepresentability
condition_ (see _e.g._ [Hastie et al., 2015]) – Lasso corresponds to having a
quadratic loss instead of $\mathcal{L}(A,\hat{\pi})$, thus
$\nabla^{2}_{A}W(A)$ in the definition below reduces to a matrix. In what
follows, we denote $u_{I}\triangleq(u_{i})_{i\in I}$ and
$U_{I,J}\triangleq(U_{i,j})_{i\in I,j\in J}$ the restriction operators.

######  Definition 1.

The _certificate_ with respect to $A$ and support $I=\\{i:A_{i}\neq 0\\}$ is

$\displaystyle
z^{\ast}_{A}\triangleq\nabla^{2}W(A)_{(:,I)}(\nabla^{2}W(A)_{(I,I)})^{-1}\operatorname{sign}(A)_{I}.$
(C)

We say that it is non-degenerate if $\|(z^{\ast}_{A})_{I^{c}}\|_{\infty}<1$.

The next proposition, whose proof can be found in Appendix A.1, shows that the
function $W(A)$ is twice differentiable, thus ensuring that C is well defined.

######  Proposition 2.

$A\mapsto W(A)$ is twice differentiable, strictly convex, with gradient and
Hessian

$\nabla_{A}W(A)=\Phi^{\ast}\pi_{A},\quad\nabla_{A}^{2}W(A)=\Phi^{\ast}\frac{\partial\pi_{A}}{\partial
A}(x,y)$

where $\pi_{A}$ is the unique solution to (1) with cost $c_{A}=\Phi A$.

The next theorem, a well-known result (see _e.g._ [Lee et al., 2015] for a
proof), shows that support recovery can be characterized via
$z^{*}_{{\widehat{A}}}$.

######  Theorem 3.

If $z^{\ast}_{{\widehat{A}}}$ is non-degenerate, then for all $\lambda$
sufficiently small, the solution $A^{\lambda}$ to
($\mathrm{iOT-}\ell_{1}(\hat{\pi})$) is sparsistent with
$\left\|A^{\lambda}-{\widehat{A}}\right\|=\mathcal{O}(\lambda)$ and
$\operatorname{Supp}\left(A^{\lambda}\right)=\operatorname{Supp}({\widehat{A}})$.

######  Remark 4.

We remark that strict convexity implies that any solution of
($\mathrm{iOT-}\ell_{1}(\hat{\pi})$) must be unique. It is then an immediate
consequence of Gamma convergence that solutions $A^{\lambda}$ to
($\mathrm{iOT-}\ell_{1}(\hat{\pi})$) converge to $\hat{A}$ as $\lambda\to 0$.
The above result shows that the support (and hence signs) are also preserved
when $\lambda$ is sufficiently small.

###  3.2 Intuition behind the certificate (C)

In this section, we provide some insight into the link between
$z^{\ast}_{{\widehat{A}}}$ and the optimization problem
($\mathrm{iOT-}\ell_{1}(\hat{\pi})$). Observe that the first order optimality
condition to ($\mathrm{iOT-}\ell_{1}(\hat{\pi})$) implies that $A^{\lambda}$
is a solution if and only if

$z^{\lambda}\triangleq-\frac{1}{\lambda}\nabla
L(A^{\lambda})\in\partial\left\|A_{\lambda}\right\|_{1},$

where the subdifferential for the $\ell_{1}$ norm has the explicit form
$\partial\left\|A\right\|_{1}=\left\\{z\;;\;\left\|z\right\|_{\infty}\leqslant
1,\;\forall
i\in\operatorname{Supp}(A),z_{i}=\operatorname{sign}(A_{i})\right\\}$. It
therefore follows $z^{\lambda}$ can be seen as a _certificate_ for the support
of $A^{\lambda}$ since
$\operatorname{Supp}(A^{\lambda})\subseteq\left\\{i\;;\;\left\lvert
z^{\lambda}_{i}\right\rvert=1\right\\}$. To study the support behaviour of
$A^{\lambda}$ for small $\lambda$, it is therefore interesting to consider the
limit of $z^{\lambda}$ as $\lambda\to 0$. In fact, its limit is precisely the
subdifferential element with minimal norm:

######  Proposition 5.

Let $z^{\lambda}\triangleq-\frac{1}{\lambda}\nabla L(A^{\lambda})$ where
$A^{\lambda}$ solves ($\mathrm{iOT-}\ell_{1}(\hat{\pi})$). Then,

$\lim_{\lambda\to
0}z^{\lambda}=z^{\min}_{{\widehat{A}}}\triangleq\operatorname*{arg\,min}_{z}\left\\{\langle
z,\,\left(\nabla^{2}W({\widehat{A}}\right)^{-1}z\rangle_{F}\;;\;z\in\partial\|{\widehat{A}}\|_{1}\right\\}$
(MNC)

Moreover, if $z^{*}_{{\widehat{A}}}$ is non-degenerate, then
$z^{\min}_{{\widehat{A}}}=z^{*}_{{\widehat{A}}}$.

The certificate (C) can therefore be seen as the _limit_ optimality condition
and determines whether the support of ${\widehat{A}}$ can be stably recovered
or not. One useful aspect of (C) is that it can be readily computed in some
cases to analyse the sparsistency properties of
($\mathrm{iOT-}\ell_{1}(\hat{\pi})$). In particular, in Section 5.1, we
provide explicit expressions for $z^{\ast}_{{\widehat{A}}}$ in the Gaussian
setting.

##  4 Sample complexity bounds

Our main contribution is the following theorem which shows that (C) also
provides a certification for sparsistency under sampling noise.

######  Theorem 6.

Let $\hat{\pi}=\mathrm{Sink}(c_{\widehat{A}},\varepsilon)$. Suppose that the
certificate $z_{{\widehat{A}}}^{*}$ is non-degenerate. Then, for all
sufficiently small regularization parameters $\lambda$ and sufficiently many
number of samples $n$,

$\lambda\lesssim
1\quad\text{and}\quad\max\left({\exp(C/\lambda)m}\lambda^{-1},\sqrt{\log(2s)}\right)\lesssim\sqrt{n},$

where $C>0$ is a constant that depends on $\hat{\pi}$, with probability at
least $1-\exp(-m^{2})$, the minimizer $A_{n}$ to ($\mathcal{P}_{n}$) is
sparsistent with ${\widehat{A}}$ and satisfies
$\left\|A_{n}-{\widehat{A}}\right\|\lesssim\lambda+\sqrt{\exp(C/\lambda)}mn^{-\tfrac{1}{2}}$.

##### Main idea behind Theorem 6

We know from Theorem 3 that there is some $\lambda_{0}>0$ such that for all
$\lambda\leqslant\lambda_{0}$, the solution to ($\mathcal{K}_{\infty}$) has
the same support as $A^{*}$. To show that the finite sample problem also
recovers the support of $A^{*}$ when $n$ is sufficiently large, we fix
$\lambda\in(0,\lambda_{0}]$ and consider the setting where the observations
are iid samples from the coupling measure $\hat{\pi}$. We will derive
convergence bounds for the primal and dual solutions as the number of samples
$n$ increases. Let $(A_{\infty},f_{\infty},g_{\infty})$ minimise
($\mathcal{K}_{\infty}$). Denote

$P_{\infty}=\frac{1}{n^{2}}(p_{\infty}(x_{i},y_{j}))_{i,j\in[n]},\quad\text{where}\quad
p_{\infty}(x,y)=\exp\left(\frac{2}{\varepsilon}\left(\Phi
A_{\infty}(x,y)+f_{\infty}(x)+g_{\infty}(y)\right)\right).$

We denote $F_{\infty}=(f_{\infty}(x_{i}))_{i\in[n]}$ and
$G_{\infty}=(g_{\infty}(y_{j}))_{j\in[n]}$. Let $(A_{n},F_{n},G_{n})$ minimise
($\mathcal{K}_{n}$), so the primal solution to ($\mathcal{P}_{n}$) is
$P_{n}\triangleq\frac{1}{n^{2}}\exp\left(\frac{2}{\varepsilon}\left(\Phi_{n}A_{n}+F_{n}\oplus
G_{n}\right)\right)$. The certificates are

$z_{\infty}=\frac{1}{\lambda}\Phi^{*}\left(p_{\infty}\alpha\otimes\beta-\hat{\pi}\right)\quad\text{and}\quad
z_{n}=\frac{1}{\lambda}\Phi_{n}^{*}\left(P_{n}-\hat{P}_{n}\right).$

By exploiting strong convexity properties of $\mathcal{J}_{n}$, one can show
the following sample complexity bound on the convergence of $z_{n}$ to
$z_{\infty}$ (the proof can be found in the appendix):

######  Proposition 7.

Let $n\gtrsim\max(m^{2}\lambda^{-2},\log(2s))$ for some $m>0$. For some
constant $C>0$, with probability at least $1-\exp(-m^{2})$

$\left\|z_{\infty}-z_{n}\right\|_{\infty}\lesssim{\exp(C/\lambda)m}{\lambda^{-1}n^{-\frac{1}{2}}}\quad\text{and}$
$\left\|A_{n}-A_{\infty}\right\|^{2}+\frac{1}{n}\sum_{i}(F_{n}-F_{\infty})_{i}^{2}+\frac{1}{n}\sum_{j}(G_{n}-G_{\infty})_{j}^{2}\lesssim\varepsilon^{2}\exp(C/\lambda)m^{2}n^{-1}.$

From Proposition 5, for all $\lambda\leqslant\lambda_{0}$ sufficiently small,
$z_{\infty}$ is nondegenerate. Moreover, the convergence result in Proposition
7 above implies that for $n$ is sufficiently large, $z_{n}$ is also
nondegenerate. Hence, since the set $\left\\{i\;;\;(z_{n})_{i}=\pm 1\right\\}$
determines the support of $A_{n}$, we have sparsistency.

##  5 Gaussian distributions

To get further insight about the sparsistency property ot iOT, we consider the
special case where the source and target distributions are Gaussians, and the
cost parametrization $c_{A}(x,y)=x^{\top}Ay$. To this end, we first derive
closed form expressions for the Hessian
$\partial_{A}^{2}\mathcal{L}(A)=\nabla_{A}^{2}W(A)$. Given
$\alpha=\mathcal{N}(m_{\alpha},\Sigma_{\alpha})$ and
$\beta=\mathcal{N}(m_{\beta},\Sigma_{\beta})$, it is known (see [Bojilov and
Galichon, 2016]) that the coupling density is also a Gaussian of the form
$\pi=\mathcal{N}\left(\binom{m_{\alpha}}{m_{\beta}},\begin{pmatrix}\Sigma_{\alpha}&\Sigma\\\
\Sigma^{\top}&\Sigma_{\beta}\end{pmatrix}\right)$ for some
$\Sigma\in\mathbb{R}^{d_{1}\times d_{2}}$. In this case, $W$ can be written as
an optimization problem over the cross covariance $\Sigma$ [Bojilov and
Galichon, 2016].

$W(A)=\sup_{\Sigma\in\mathbb{R}^{d_{1}\times d_{2}}}\langle
A,\,\Sigma\rangle+\frac{\varepsilon}{2}\log\det\left(\Sigma_{\beta}-\Sigma^{\top}\Sigma_{\alpha}^{-1}\Sigma\right),$
(3)

with the optimal solution being precisely the cross-covariance of the optimal
coupling $\pi$. In [Bojilov and Galichon, 2016], the authors provide an
explicit formula for the minimizer $\Sigma$ to (3), and, consequently, for
$\nabla W(A)$:

$\Sigma=\Sigma_{\alpha}A\Delta\left(\Delta
A^{\top}\Sigma_{\alpha}A\Delta\right)^{-\frac{1}{2}}\Delta-\tfrac{1}{2}\varepsilon
A^{\dagger,\top}\quad\text{where}\quad\Delta\triangleq\left(\Sigma_{\beta}+\frac{\varepsilon^{2}}{4}A^{\dagger}\Sigma_{\alpha}^{-1}A^{\dagger,\top}\right)^{\frac{1}{2}}.$
(4)

By differentiating the first order condition for $W$, that is
$A^{\dagger}=\varepsilon^{-1}(\Sigma_{\beta}-\Sigma^{\top}\Sigma_{\alpha}^{-1}\Sigma)\Sigma^{\dagger}\Sigma_{\alpha},$
Galichon derives the expression for the Hessian in terms of $\Sigma$ in (4):

$\nabla^{2}W(A)=\partial_{A}\Sigma=\varepsilon\left(\Sigma_{\alpha}\Sigma^{-1,\top}\otimes\Sigma_{\beta}\Sigma^{-1}+\mathbb{T}\right)^{-1}\left(A^{-1,\top}\otimes
A^{-1}\right),$ (5)

where $\mathbb{T}$ is such that
$\mathbb{T}\mathrm{vec}(A)=\mathrm{vec}(A^{\top})$. This formula does not hold
in when $A$ is rectangular or rank deficient, since $A^{\dagger}$ is not
differentiable. In the following, we derive, via the implicit function
theorem, a general formula for $\partial_{A}\Sigma$ that agrees with that of
Galichon in the square invertible case.

######  Lemma 8.

Denoting $\Sigma$ as in (4), one has

$\displaystyle\nabla^{2}W(A)$
$\displaystyle=\varepsilon\left(\varepsilon^{2}(\Sigma_{\beta}-\Sigma^{\top}\Sigma_{\alpha}\Sigma)^{-1}\otimes(\Sigma_{\alpha}-\Sigma\Sigma_{\beta}^{-1}\Sigma^{\top})^{-1}+(A^{\top}\otimes
A)\mathbb{T}\right)^{-1}$ (6)

This formula given in Lemma 8 provides an explicit expression for the
certificate (C).

###  5.1 Limit cases for large and small $\varepsilon$

This section explores the behaviour of the certificate in the large/small
$\varepsilon$ limits: Proposition 9 reveals that the large $\varepsilon$ limit
coincides with the classical Lasso while Proposition 10 reveals that the small
epsilon limit (for symmetric $A\succ 0$ and
$\Sigma_{\alpha}=\Sigma_{\beta}=\mathrm{Id}$) coincides with the Graphical
Lasso. In the following results, we denote the functional in
($\mathrm{iOT-}\ell_{1}(\hat{\pi})$) with parameters $\lambda$ and
$\varepsilon$ by $\mathcal{F}(A)_{\varepsilon,\lambda}$.

######  Proposition 9 ($\varepsilon\to\infty$).

Let ${\widehat{A}}$ be invertible and let
$\hat{\pi}=\mathrm{Sink}(c_{\widehat{A}},\varepsilon)$ be the observed
coupling between $\alpha=\mathcal{N}(m_{\alpha},\Sigma_{\alpha})$ and
$\beta=\mathcal{N}(m_{\beta},\Sigma_{\beta})$. Then,

$\lim_{\varepsilon\to\infty}z_{\varepsilon}=(\Sigma_{\beta}\otimes\Sigma_{\alpha})_{(:,I)}\big{(}(\Sigma_{\beta}\otimes\Sigma_{\alpha})_{(I,I)}\big{)}^{-1}\operatorname{sign}({\widehat{A}})_{I}.$
(7)

Moreover, for $\lambda_{0}>0$, given any sequence $(\varepsilon_{j})_{j}$ and
$A_{j}\in\operatorname*{arg\,min}_{A}\mathcal{F}(A)_{\varepsilon,\lambda_{0}/\varepsilon_{j}}(A)$
with $\lim_{j\to\infty}\varepsilon_{j}=\infty$, any cluster point of
$(A_{j})_{j}$ is in

$\underset{A\in\mathbb{R}^{d\times
d}}{\operatorname*{arg\,min}}\;\lambda_{0}\left\|A\right\|_{1}+\frac{1}{2}\|{(\Sigma_{\beta}^{1/2}\otimes\Sigma_{\alpha}^{1/2})(A-{\widehat{A}})}\|_{F}^{2}$
(8)

######  Proposition 10 ($\varepsilon\to 0$).

Let ${\widehat{A}}$ be symmetric positive-definite and let
$\hat{\pi}=\mathrm{Sink}(c_{\widehat{A}},\varepsilon)$ be the observed
coupling between $\alpha=\mathcal{N}(m_{\alpha},\mathrm{Id})$ and
$\beta=\mathcal{N}(m_{\beta},\mathrm{Id})$. Then,

$\lim_{\varepsilon\to
0}z_{\varepsilon}=({\widehat{A}}^{-1}\otimes{\widehat{A}}^{-1})_{(:,I)}\big{(}({\widehat{A}}^{-1}\otimes{\widehat{A}}^{-1})_{(I,I)}\big{)}^{-1}\operatorname{sign}({\widehat{A}})_{I}.$
(9)

Let $\lambda_{0}>0$. Then, optimizing over symmetric positive semi-definite
matrices, given any sequence $(\varepsilon_{j})_{j}$ and
$A_{j}\in\operatorname*{arg\,min}_{A\succeq
0}\mathcal{F}(A)_{\varepsilon,\lambda_{0}\varepsilon_{j}}(A)$ with
$\lim_{j\to\infty}\varepsilon_{j}=0$, any cluster point of $(A_{j})_{j}$ is in

$\underset{A\succeq
0}{\operatorname*{arg\,min}}\;\lambda_{0}\left\|A\right\|_{1}-\frac{1}{2}\log\det(A)+\frac{1}{2}\langle
A,\,{\widehat{A}}^{-1}\rangle.$ (10)

Observe that the $\varepsilon\to\infty$ and $\varepsilon\to 0$ limits given in
(7) and (9) correspond respectively to the certificates associated to a Lasso
optimization problem, _i.e.,_ problem (8), and to a graphical Lasso
optimization problem, _i.e.,_ problem 10. This can be easily derived by
plugging in the Hessian of the unregularized loss functions of (8) and (10)
into (C); alternatively, see [Hastie et al., 2015] for the Lasso certificate
and [Ravikumar et al., 2011] for graphical Lasso.

###  5.2 Numerical illustrations

In order to gain some insight into the impact of $\varepsilon$ and the
covariance structure on the efficiency of iOT, we present numerical
computations of certificates here. We fix the covariances of the input
measures as $\Sigma_{\alpha}=\Sigma_{\beta}=\mathrm{Id}_{n}$, 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=\delta\mathrm{Id}_{n}+\operatorname{diag}(G1_{n})-G\in\mathbb{R}^{n\times
n}$ is defined as a shifted Laplacian matrix of some graph adjacency matrix
$G$, for a graph of size $n=80$ (similar conclusions hold for larger graphs).
We set the shift $\delta$ 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 $\varepsilon$, 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
a spatially smoothly varying covariance $C^{-1}$. Figure 1 illustrates how the
value of the certificates $z_{i,j}$ evolves depending on the indexes $(i,j)$
for three types of graphs (circular, planar, and Erdős–Rényi with a
probability of edges equal to 0.1), for several values of $\varepsilon$. By
construction, $z_{i,i}=1$ and $z_{i,j}=-1$ for $(i,j)$ connected by the graph.
For $z$ to be non-degenerate, it is required that $|z_{i,j}|<1$, as $i$ moves
away from $j$ on the graph. For the circular and planar graphs, the horizontal
axis represents the geodesic distance $d_{\text{geod}}(i,j)$, demonstrating
how the certificates become well-behaved as the distance increases. The planar
graph displays envelope curves showing the range of values of $z$ for a fixed
value of $d_{\text{geod}}(i,j)$, while this is a single-valued curve for the
circular graph due to periodicity. For the Erdős–Rényi graph, to better
account for the randomness of the certificates along the graph edges, we
display the histogram of the distribution of $|z_{i,j}|$ for
$d_{\text{geod}}(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 $\varepsilon$ improves the behavior of the
certificates (which is in line with the stability analysis of Section 5.1),
and that pairs of vertices $(i,j)$ connected by a small distance
$d_{\text{geod}}(i,j)$ are the most likely to be degenerate. This suggests
that they will be inaccurately estimated by iOT for small $\varepsilon$.

|  |
---|---|---
|  |
Circular | Planar | Erdős–Rényi
Figure 1: Display of the certificate values $z_{i,j}$ for three types of
graphs, for varying $\varepsilon$. Left, middle: plotted as a function of the
geodesic distance $d_{\text{geod}}(i,j)$ on the $x$-axis. Right: histogram of
$z_{i,j}$ for $(i,j)$ at distance $d_{\text{geod}}(i,j)=2$.

Figure 2 displays the recovery performances of $\ell^{1}$–iOT for the circular
graph shown on the left of Figure 1 (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 D. The
performance is represented using the number of inaccurately estimated
coordinates in the estimated cost matrix $C$, so a score of 0 means a perfect
estimation of the support (sparsistency is achieved). For $\varepsilon=0.1$,
sparsistency cannot be achieved, aligning with the fact that the certificate
$z$ is degenerate as depicted in the previous figure. In sharp contrast, for
larger $\varepsilon$, sparsistency can be attained as soon as the number of
samples $N$ is sufficiently large, which also aligns with our theory of
sparsistency of $\ell^{1}$–iOT.

|  |
---|---|---
$\varepsilon=0.1$ | $\varepsilon=1$ | $\varepsilon=10$
Figure 2:  Recovery performance (number of wrongly estimated position) of
$\ell^{1}$–iOT as a function of $\lambda$ for three different values of
$\varepsilon$.

Figure 3 displays the certificate values in the case of a non-symmetric planar
graph. The graph is obtained by deleting all edges on a planar graph with
$i<j$. We plot the certificate values as a function of geodesic distances
$d_{\text{geod}}(i,j)$ of the symmetrized graph. The middle plot shows the
certificate values on $i\geqslant j$ (where the actual edges are constrained
to be and nondegeneracy requires values smaller than 1 in absolute value for
$d_{\text{geod}}(i,j)\geqslant 2$). The right plot shows the certificate
values on $i\leqslant j$ where there are no edges and for nondegeneracy, one
expects values smaller than 1 in absolute value for
$d_{\text{geod}}(i,j)\geqslant 1$. Observe that here, the certificate is
degenerate for small values of $d_{\text{geod}}(i,j)$ when $\varepsilon=0$,
meaning that the problem is unstable at the “ghost” symmetric edges. As
$\varepsilon\to\infty$, the certificate becomes non-degenerate. Observe that
as $\varepsilon\to\infty$, the certificate becomes non-degenerate.

|
---|---
Figure 3:  Display of certificate values for a non-symmetric planar graph, for
varying $\varepsilon$ with edges only for $i>j$. Middle/Right: plots of the
certificate values as a function of the _geodesic distance
$d_{\text{geod}}(i,j)$ of the symmetrized graph_. The middle plot show the
values when restricted to $i\geqslant j$. The right plot shows the values
restricted to $i\leqslant j$ (where there are no edges).

## Conclusion

In this paper, we have proposed the first theoretical analysis of the recovery
performance of $\ell^{1}$-iOT. Much of 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-posed and fails to maintain sparsity for
overly small $\varepsilon$. 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’avenir” program, reference
ANR19-P3IA-0001 (PRAIRIE 3IA Institute).

## References

  * [Ahuja and Orlin, 2001] Ahuja, R. K. and Orlin, J. B. (2001).  Inverse optimization.  Operations research, 49(5):771–783.
  * [Altschuler et al., 2017] Altschuler, J., Niles-Weed, J., and Rigollet, P. (2017).  Near-linear time approximation algorithms for optimal transport via sinkhorn iteration.  Advances in neural information processing systems, 30.
  * [Arjovsky et al., 2017] Arjovsky, M., Chintala, S., and Bottou, L. (2017).  Wasserstein generative adversarial networks.  In International conference on machine learning, pages 214–223. PMLR.
  * [Bellet et al., 2013] Bellet, A., Habrard, A., and Sebban, M. (2013).  A survey on metric learning for feature vectors and structured data.  arXiv preprint arXiv:1306.6709.
  * [Bojilov and Galichon, 2016] Bojilov, R. and Galichon, A. (2016).  Matching in closed-form: equilibrium, identification, and comparative statics.  Economic Theory, 61:587–609.
  * [Borwein and Lewis, 2006] Borwein, J. and Lewis, A. (2006).  Convex Analysis.  Springer.
  * [Candès et al., 2006] Candès, E. J., Romberg, J., and Tao, T. (2006).  Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information.  IEEE Transactions on information theory, 52(2):489–509.
  * [Carlier et al., 2023] Carlier, G., Dupuy, A., Galichon, A., and Sun, Y. (2023).  Sista: learning optimal transport costs under sparsity constraints.  Communications on Pure and Applied Mathematics, 76(9):1659–1677.
  * [Chiu et al., 2022] Chiu, W.-T., Wang, P., and Shafto, P. (2022).  Discrete probabilistic inverse optimal transport.  In International Conference on Machine Learning, pages 3925–3946. PMLR.
  * [Colson et al., 2005] Colson, B., Marcotte, P., and Savard, G. (2005).  Bilevel programming: A survey.  4or, 3:87–107.
  * [Cuturi, 2013] Cuturi, M. (2013).  Sinkhorn distances: Lightspeed computation of optimal transport.  In Adv. in Neural Information Processing Systems, pages 2292–2300.
  * [Cuturi and Avis, 2014] Cuturi, M. and Avis, D. (2014).  Ground metric learning.  The Journal of Machine Learning Research, 15(1):533–564.
  * [Cuturi and Peyré, 2016] Cuturi, M. and Peyré, G. (2016).  A smoothed dual approach for variational wasserstein problems.  SIAM Journal on Imaging Sciences, 9(1):320–343.
  * [Davis and Dhillon, 2008] Davis, J. V. and Dhillon, I. S. (2008).  Structured metric learning for high dimensional problems.  In Proceedings of the 14th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 195–203.
  * [De Bortoli et al., 2021] De Bortoli, V., Thornton, J., Heng, J., and Doucet, A. (2021).  Diffusion schrödinger bridge with applications to score-based generative modeling.  Advances in Neural Information Processing Systems, 34:17695–17709.
  * [Deimling, 2010] Deimling, K. (2010).  Nonlinear functional analysis.  Courier Corporation.
  * [Dupuy and Galichon, 2014] Dupuy, A. and Galichon, A. (2014).  Personality traits and the marriage market.  Journal of Political Economy, 122(6):1271–1319.
  * [Dupuy et al., 2019] Dupuy, A., Galichon, A., and Sun, Y. (2019).  Estimating matching affinity matrices under low-rank constraints.  Information and Inference: A Journal of the IMA, 8(4):677–689.
  * [Eisenberger et al., 2022] Eisenberger, M., Toker, A., Leal-Taixé, L., Bernard, F., and Cremers, D. (2022).  A unified framework for implicit sinkhorn differentiation.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 509–518.
  * [Galichon, 2018] Galichon, A. (2018).  Optimal transport methods in economics.  Princeton University Press.
  * [Galichon and Salanié, 2010] Galichon, A. and Salanié, B. (2010).  Matching with trade-offs: Revealed preferences over competing characteristics.  Preprint hal-00473173.
  * [Galichon and Salanié, 2022] Galichon, A. and Salanié, B. (2022).  Cupid’s invisible hand: Social surplus and identification in matching models.  The Review of Economic Studies, 89(5):2600–2629.
  * [Genevay et al., 2019] Genevay, A., Chizat, L., Bach, F., Cuturi, M., and Peyré, G. (2019).  Sample complexity of Sinkhorn divergences.  In The 22nd international conference on artificial intelligence and statistics, pages 1574–1583. PMLR.
  * [Genevay et al., 2018] Genevay, A., Peyré, G., and Cuturi, M. (2018).  Learning generative models with sinkhorn divergences.  In International Conference on Artificial Intelligence and Statistics, pages 1608–1617. PMLR.
  * [Hastie et al., 2015] Hastie, T., Tibshirani, R., and Wainwright, M. (2015).  Statistical learning with sparsity: the lasso and generalizations.  CRC press.
  * [Heitz et al., 2021] Heitz, M., Bonneel, N., Coeurjolly, D., Cuturi, M., and Peyré, G. (2021).  Ground metric learning on graphs.  Journal of Mathematical Imaging and Vision, 63:89–107.
  * [Hiriart-Urruty et al., 1993] Hiriart-Urruty, J.-B., Lemaréchal, C., Hiriart-Urruty, J.-B., and Lemaréchal, C. (1993).  Conjugacy in convex analysis.  Convex Analysis and Minimization Algorithms II: Advanced Theory and Bundle Methods, pages 35–90.
  * [Huang et al., 2016] Huang, G., Guo, C., Kusner, M. J., Sun, Y., Sha, F., and Weinberger, K. Q. (2016).  Supervised word mover’s distance.  Advances in neural information processing systems, 29.
  * [Kerdoncuff et al., 2021] Kerdoncuff, T., Emonet, R., and Sebban, M. (2021).  Metric Learning in Optimal Transport for Domain Adaptation.  In International Joint Conference on Artificial Intelligence, Kyoto, Japan.
  * [Kusner et al., 2015] Kusner, M., Sun, Y., Kolkin, N., and Weinberger, K. Q. (2015).  From word embeddings to document distances.  In Proc. of the 32nd Intern. Conf. on Machine Learning, pages 957–966.
  * [Lee et al., 2015] Lee, J. D., Sun, Y., and Taylor, J. E. (2015).  On model selection consistency of regularized m-estimators.  Electronic Journal of Statistics, 9:608–642.
  * [Léonard, 2012] Léonard, C. (2012).  From the Schrödinger problem to the Monge–Kantorovich problem.  Journal of Functional Analysis, 262(4):1879–1920.
  * [Li et al., 2019] Li, R., Ye, X., Zhou, H., and Zha, H. (2019).  Learning to match via inverse optimal transport.  Journal of machine learning research, 20.
  * [Ma et al., 2020] Ma, S., Sun, H., Ye, X., Zha, H., and Zhou, H. (2020).  Learning cost functions for optimal transport.  arXiv preprint arXiv:2002.09650.
  * [Nutz, 2021] Nutz, M. (2021).  Introduction to entropic optimal transport.  Lecture notes, Columbia University.
  * [Paty and Cuturi, 2020] Paty, F.-P. and Cuturi, M. (2020).  Regularized optimal transport is ground cost adversarial.  In International Conference on Machine Learning, pages 7532–7542. PMLR.
  * [Peyré et al., 2019] Peyré, G., Cuturi, M., et al. (2019).  Computational optimal transport: With applications to data science.  Foundations and Trends® in Machine Learning, 11(5-6):355–607.
  * [Poon and Peyré, 2021] Poon, C. and Peyré, G. (2021).  Smooth bilevel programming for sparse regularization.  Advances in Neural Information Processing Systems, 34:1543–1555.
  * [Ravikumar et al., 2011] Ravikumar, P., Wainwright, M. J., Raskutti, G., and Yu, B. (2011).  High-dimensional covariance estimation by minimizing l1-penalized log-determinant divergence.  Electronic Journal of Statistics, 5:935–980.
  * [Rigollet and Stromme, 2022] Rigollet, P. and Stromme, A. J. (2022).  On the sample complexity of entropic optimal transport.  arXiv preprint arXiv:2206.13472.
  * [Rubner et al., 2000] Rubner, Y., Tomasi, C., and Guibas, L. J. (2000).  The earth mover’s distance as a metric for image retrieval.  International Journal of Computer Vision, 40(2):99–121.
  * [Santambrogio, 2015] Santambrogio, F. (2015).  Optimal transport for applied mathematicians.  Birkäuser, NY, 55(58-63):94.
  * [Schiebinger et al., 2019] Schiebinger, G., Shu, J., Tabaka, M., Cleary, B., Subramanian, V., Solomon, A., Gould, J., Liu, S., Lin, S., Berube, P., et al. (2019).  Optimal-transport analysis of single-cell gene expression identifies developmental trajectories in reprogramming.  Cell, 176(4):928–943.
  * [Shi et al., 2023] Shi, L., Zhang, G., Zhen, H., Fan, J., and Yan, J. (2023).  Understanding and generalizing contrastive learning from the inverse optimal transport perspective.
  * [Sinkhorn, 1964] Sinkhorn, R. (1964).  A relationship between arbitrary positive matrices and doubly stochastic matrices.  Ann. Math. Statist., 35:876–879.
  * [Stuart and Wolfram, 2020] Stuart, A. M. and Wolfram, M.-T. (2020).  Inverse optimal transport.  SIAM Journal on Applied Mathematics, 80(1):599–619.
  * [Tibshirani, 1996] Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 58(1):267–288.
  * [Tropp et al., 2015] Tropp, J. A. et al. (2015).  An introduction to matrix concentration inequalities.  Foundations and Trends® in Machine Learning, 8(1-2):1–230.
  * [Vaiter et al., 2015] Vaiter, S., Golbabaee, M., Fadili, J., and Peyré, G. (2015).  Model selection with low complexity priors.  Information and Inference: A Journal of the IMA, 4(3):230–287.
  * [Wang and Guibas, 2012] Wang, F. and Guibas, L. J. (2012).  Supervised earth mover’s distance learning and its computer vision applications.  In European Conference on Computer Vision, pages 442–455. Springer.
  * [Weinberger et al., 2006] Weinberger, K. Q., Blitzer, J., and Saul, L. K. (2006).  Distance metric learning for large margin nearest neighbor classification.  In Advances in neural information processing systems, pages 1473–1480.
  * [Xing et al., 2002] Xing, E., Jordan, M., Russell, S. J., and Ng, A. (2002).  Distance metric learning with application to clustering with side-information.  Advances in neural information processing systems, 15.
  * [Xu et al., 2018] Xu, J., Luo, L., Deng, C., and Huang, H. (2018).  Multi-level metric learning via smoothed wasserstein distance.  In IJCAI, pages 2919–2925.
  * [Zen et al., 2014] Zen, G., Ricci, E., and Sebe, N. (2014).  Simultaneous ground metric learning and matrix factorization with earth mover’s distance.  In 2014 22nd International Conference on Pattern Recognition, pages 3690–3695. IEEE.
  * [Zhang and Liu, 1996] Zhang, J. and Liu, Z. (1996).  Calculating some inverse linear programming problems.  Journal of Computational and Applied Mathematics, 72(2):261–273.
  * [Zhao and Yu, 2006] Zhao, P. and Yu, B. (2006).  On model selection consistency of lasso.  The Journal of Machine Learning Research, 7:2541–2563.

##  Appendix A Proofs for Section 3

###  A.1 Proof of Proposition 2 (differentiability of $W$)

For the strict convexity of $W(A)$ see Lemma 3 in [Dupuy and Galichon, 2014].
The gradient formula can also be found in [Dupuy and Galichon, 2014] and
trivially follows from the envelope theorem because the optimization problem
$W(A)$ has a unique solution. We will only give a proof of the Hessian
formula.

The formula for the Hessian follows from the formula for the gradient provided
we show that the density $\pi_{A}$ is continuously differentiable with respect
to $A$. The fact that we can swap the order of the operator $\Phi^{\ast}$ and
partial differentiation follows from the conditions for differentiability
under the integral sign which hold since the measures are compactly supported.
Without loss of generality, let $\varepsilon=1$. Then, the optimizer in $W(A)$
is of the form [Santambrogio, 2015]:

$\displaystyle\pi_{A}(x,y)=\exp\Big{(}u_{A}(x)+v_{A}(y)+c_{A}(x,y)\Big{)},$

where $u_{A}(\cdot)$ and $v_{A}(\cdot)$ satisfy

$\displaystyle u_{A}(x)$
$\displaystyle=-\log\int_{\mathcal{Y}}\exp\Big{(}v_{A}(y)+c_{A}(x,y)\Big{)}\,d\beta(y),\quad\alpha\text{-a.s.}$
$\displaystyle v_{A}(y)$
$\displaystyle=-\log\int_{\mathcal{X}}\exp\Big{(}u_{A}(x)+c_{A}(x,y)\Big{)}\,d\alpha(x),\quad\beta\text{-a.s.}$

It is known (see _e.g._ [Nutz, 2021]) that the functions $u_{A}(\cdot)$ and
$v_{A}(\cdot)$ inherit the modulus of continuity of the cost (in this case the
map $(x,y)\to c_{A}(x,y))$) and, hence, since we are assuming the measures to
be compactly supported, it follows that $u_{A}\in L^{2}(\alpha)$ and $v_{A}\in
L^{2}(\beta)$. Moreover, if $(u_{A},v_{A})$ solve these two equations then,
for any constant $c$, the pair $(u_{A}+c,v_{A}-c)$ is also a solution and,
hence, to eliminate this ambiguity we consider solutions in
$\mathcal{H}=L^{2}(\alpha)\times L_{0}^{2}(\beta)$, where
$L_{0}^{2}(\beta)=\left\\{g\in L^{2}(\beta)\;;\;\int g\,d\beta(y)=0\right\\}$.

To show that $W$ is twice differentiable, since $\nabla W(A)=\Phi^{*}\pi_{A}$,
it is sufficient to show that $A\mapsto u_{A}$ and $A\mapsto v_{A}$ are
differentiable. To this end, we will apply the Implicit Function Theorem (in
Banach spaces) to the map

$\displaystyle F:\mathcal{H}\times\mathbb{R}^{s}$ $\displaystyle\to
L^{2}(\alpha)\times L^{2}(\beta)$ $\displaystyle\begin{bmatrix}u\\\ v\\\
A\end{bmatrix}$
$\displaystyle\to\begin{bmatrix}u+\log\int_{\mathcal{Y}}\exp\Big{(}v(y)+c_{A}(x,y)\Big{)}\,d\beta(y)\\\
v+\log\int_{\mathcal{X}}\exp\Big{(}u(x)+c_{A}(x,y)\Big{)}\,d\alpha(x)\end{bmatrix},$

since we have $F(u_{A},v_{A},A)=0$. The partial derivative of $F$ at
$(u_{A},v_{A},A)$ , denoted by $\partial_{u,v}F(u_{A},v_{A},A)$, is the linear
map defined by

$\displaystyle\Big{(}\partial_{u,v}F(u_{A},v_{A},A)\Big{)}(f,g)=\Big{(}f+\int_{\mathcal{Y}}p_{A}(\cdot,y)g(y)\,d\beta(y),g+\int_{\mathcal{X}}q_{A}(x,\cdot)f(x)\,d\alpha(x)\Big{)},$
(11)

where

$\displaystyle
p_{A}(x,y)=\frac{\exp\Big{(}v_{A}(y)+c_{A}(x,y)\Big{)}}{\int_{\mathcal{Y}}\exp\Big{(}v_{A}(y)+c_{A}(x,y)\Big{)}\,d\beta(y)}\text{
and
}q_{A}(x,y)=\frac{\exp\Big{(}u_{A}(x)+c_{A}(x,y)\Big{)}}{\int_{\mathcal{X}}\exp\Big{(}u_{A}(y)+c_{A}(x,y)\Big{)}\,d\beta(x)}.$

Note that, since $F(u_{A},v_{A},A)=0$, $p_{A}(x,y)=q_{A}(x,y)=\pi_{A}(x,y)$.
Moreover, since $\pi_{A}$ has marginals $\alpha$ and $\beta$, it follows that

$\begin{aligned}
&\Big{\langle}(f,g),\Big{(}\partial_{u,v}F(u_{A},v_{A},A)\Big{)}(f,g)\Big{\rangle}\\\
=&\int_{\mathcal{X}}f(x)^{2}\,d\alpha(x)+\int_{\mathcal{X}\times\mathcal{Y}}2f(x)g(y)\pi_{A}(x,y)\,d(\alpha\otimes\beta)(x,y)+\int_{\mathcal{Y}}g(y)^{2}\,d\beta(y)\\\
=&\int_{\mathcal{X}\times\mathcal{Y}}\Big{(}f(x)+g(y)\Big{)}^{2}\pi_{A}(x,y)\,d(\alpha\otimes\beta)(x,y)\end{aligned}.$
(12)

This shows that $\partial_{u,v}F(u_{A},v_{A},A)$ is invertible – the last is
line of LABEL:InnerProductPartialF is zero if and only if $f\oplus g\equiv 0$
and since $g\in L^{2}_{0}(\beta)$ it follows that $g=0$ and $f=0$.

To conclude the proof we need to show that
$\left(\partial_{u,v}F(u_{A},v_{A},A)\right)^{-1}$ is a bounded operator (see
_e.g._ [Deimling, 2010] for the statement of IFT in Banach spaces) and to show
this it is enough to show that, for some constant $C$,

$\Big{\|}\Big{(}\partial_{u,v}F(u_{A},v_{A},A)\Big{)}(f,g)\Big{\|}\geqslant
C\|(f,g)\|.$

This follows from LABEL:InnerProductPartialF and the fact that there exists a
constant $C$ such that $\pi_{A}(x,y)\geqslant C$ for all $x$ and $y$ (see
_e.g._ [Nutz, 2021] for the existence of $C$). In fact, from $g\in
L^{2}_{0}(\beta)$, we obtain
$\int(f(x)+g(y)))^{2}\,d(\alpha\otimes\beta)(x,y)=\|(f,g)\|^{2}$ and, hence,
LABEL:InnerProductPartialF implies that

$\displaystyle\Big{\langle}(f,g),\Big{(}\partial_{u,v}F(u_{A},v_{A},A)\Big{)}(f,g)\Big{\rangle}\geqslant
C\|(f,g)\|^{2}$

and from Cauchy Schwartz applied to the left-hand-side we obtain

$\displaystyle\Big{\|}\Big{(}\partial_{u,v}F(u_{A},v_{A},A)\Big{)}(f,g)\Big{\|}\geqslant
C\|(f,g)\|,$

thus concluding the proof.

###  A.2 Connection of the certificate (C) with optimization problem
($\mathrm{iOT-}\ell_{1}(\hat{\pi})$)

As mentioned in Section 3.2, the vector $z^{\lambda}$ provides insight into
support recovery since
$\operatorname{Supp}(A^{\lambda})\subseteq\left\\{i\;;\;z^{\lambda}_{i}=\pm
1\right\\}$. In this section, we provide the proof to Proposition 5, which
shows that as $\lambda$ converges to $0$, $z^{\lambda}$ converges to the
solution of a quadratic optimization problem (MNC) that we term the _minimal
norm certificate_.

Note that the connection with (C) 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 ${\widehat{A}}$) in
(MNC) can be dropped since they are inactive; in this case (MNC) reduces to a
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 is clear. So, under non-degeneracy, (C) can be seen as the _limit_
optimality vector and determines the support of $A^{\lambda}$ when $\lambda$
is small.

To prove Proposition 5, we first show that the vector $z^{\lambda}$ coincides
the solution to a dual problem of ($\mathrm{iOT-}\ell_{1}(\hat{\pi})$).

######  Proposition 11.

Let $W^{*}$ be the convex conjugate of $W$. Problem
($\mathrm{iOT-}\ell_{1}(\hat{\pi})$) admits a dual given by

$\operatorname*{arg\,min}_{z}W^{\ast}(\hat{\Sigma}_{xy}-\lambda z)\quad\text{
subject to }\quad\|z\|_{\infty}\leqslant 1,$ (13)

where $\hat{\Sigma}_{xy}=\Phi^{\ast}\hat{\pi}$ . Moreover, a pair of primal-
dual solutions $(A^{\lambda},z^{\lambda})$ is related by

$z^{\lambda}=-\tfrac{1}{\lambda}\nabla_{A}\mathcal{L}(A^{\lambda},\hat{\pi})\quad\text{and}\quad
z^{\lambda}\in\partial\|{A^{\lambda}}\|_{1}.$ (14)

###### Proof.

Observe that we can write $\mathcal{L}(A,\hat{\pi})$ as

$\displaystyle\mathcal{L}(A,\hat{\pi})$
$\displaystyle=W(A)-\int_{\mathcal{X}\times\mathcal{Y}}(\Phi
A)(x,y)\,d\hat{\pi}=W(A)-\langle\Phi^{\ast}\hat{\pi},A\rangle_{F}=W(A)-\langle
A,\hat{\Sigma}_{xy}\rangle.$

The _Fenchel Duality Theorem_ (see e.g. [Borwein and Lewis, 2006]) yields a
dual of ($\mathrm{iOT-}\ell_{1}(\hat{\pi})$) given by

$\displaystyle\operatorname*{arg\,min}_{w}W^{\ast}(\hat{\Sigma}_{xy}-w)+(\lambda\|\cdot\|_{1})^{\ast}(w).$

To conclude the proof just note that the Fenchel conjugate of
$\lambda\|\cdot\|_{1}$ is the indicator of the set
$\\{v:\|v\|_{\infty}\leqslant\lambda\\}$ and make a change of variable
$z\triangleq 1/\lambda w$ to obtain 13. The relationship between any primal-
dual pair in Fenchel Duality can also be found in [Borwein and Lewis, 2006].

∎

###### Proof of Proposition 5.

We begin by noting that $W^{\ast}$ is of class $C^{2}$ in a neighborhood of
$\Sigma_{xy}$ and that

$\displaystyle\nabla^{2}W^{\ast}(\hat{\Sigma}_{xy})=\Big{(}\nabla^{2}W({\widehat{A}})\Big{)}^{-1}.$
(15)

To see this, note that Proposition 2 together with the assumption on
$\hat{\pi}$ implies that

$\displaystyle\hat{\Sigma}_{xy}=\nabla W({\widehat{A}}).$

Moreover, since $W(A)$ is twice continuously differentiable and strictly
convex (see Proposition C), it follows (see _e.g._ [Hiriart-Urruty et al.,
1993]) that $W^{\ast}(\cdot)$ is $C^{2}$ and strictly convex in a neighborhood
of $\hat{\Sigma}_{xy}$ and that (15) holds.

Now observe that, since $\nabla W^{\ast}(\hat{\Sigma}_{xy})=\nabla
W^{\ast}\big{(}\nabla W({\widehat{A}})\big{)}={\widehat{A}}$, we can rewrite
$\partial\left\|{\widehat{A}}\right\|_{1}$ as

$\partial\left\|{\widehat{A}}\right\|_{1}=\operatorname*{arg\,min}_{z}\big{\langle}-z,\nabla
W^{\ast}(\hat{\Sigma}_{xy})\big{\rangle}\quad\text{subject
to}\quad\|z\|_{\infty}\leqslant 1$ (16)

Observe that since $z^{ambda}$ are uniform bounded vectors due to the
constraint set in 13, there is a convergent subsequence converging to some
$z^{\ast}$. We later deduce that all limit points are the same and hence, the
full sequence $z^{\lambda}$ converges to $z^{\ast}$. Let $\lambda_{n}$ be such
that $\lim_{\lambda_{n}\to 0}z^{\lambda_{n}}=z^{\ast}$, and let $z^{0}$ be any
element in $\partial\left\|{\widehat{A}}\right\|_{1}$. We have that

$\displaystyle\big{\langle}-z^{0},\nabla
W^{\ast}(\hat{\Sigma}_{xy})\big{\rangle}$
$\displaystyle\leqslant\big{\langle}-z^{\lambda_{n}},\nabla
W^{\ast}(\hat{\Sigma}_{xy})\big{\rangle}\leqslant\frac{1}{\lambda_{n}}\Big{(}W^{\ast}\big{(}\hat{\Sigma}_{xy}-\lambda_{n}z^{\lambda_{n}}\big{)}-W(\hat{\Sigma}_{xy})\Big{)}$
$\displaystyle\leqslant\frac{1}{\lambda_{n}}\Big{(}W^{\ast}\big{(}\hat{\Sigma}_{xy}-\lambda_{n}z^{0}\big{)}-W(\hat{\Sigma}_{xy})\Big{)},$

where the first inequality is the optimality of $z^{0}$, the second inequality
is the gradient inequality of convex functions and the last inequality follows
from optimality of $z^{\lambda}$. Taking the limit as $\lambda_{n}\to 0$ we
obtain that

$\displaystyle\big{\langle}-z^{0},\nabla
W^{\ast}(\hat{\Sigma}_{xy})\big{\rangle}=\big{\langle}-z^{\ast},\nabla
W^{\ast}(\hat{\Sigma}_{xy})\big{\rangle},$

showing that $z^{\ast}\in\partial\left\|{\widehat{A}}\right\|_{1}$. We now
finish the proof by showing that

$\displaystyle\Big{\langle}z^{\ast},\nabla^{2}W^{\ast}(\hat{\Sigma}_{xy})z^{\ast}\Big{\rangle}\leqslant\Big{\langle}z^{0},\nabla^{2}W^{\ast}(\hat{\Sigma}_{xy})z^{0}\Big{\rangle}.$
(17)

Since $W^{\ast}(\cdot)$ is $C^{2}$ in a neighborhood of $\hat{\Sigma}_{xy}$,
Taylor’s theorem ensures that there exists a remainder function $R(x)$ with
$\lim_{x\to\hat{\Sigma}_{xy}}R(x)=0$ such that

$\displaystyle\langle-z^{\lambda_{n}},\nabla
W^{\ast}(\hat{\Sigma}_{xy})\rangle+\frac{\lambda_{n}}{2}\big{\langle}z^{\lambda_{n}},\nabla^{2}W^{\ast}(\hat{\Sigma}_{xy})z^{\lambda_{n}}\big{\rangle}+R(\hat{\Sigma}_{xy}+\lambda_{n}z^{\lambda_{n}})\lambda^{2}_{n}$
$\displaystyle=$
$\displaystyle\frac{1}{\lambda_{n}}\Big{(}W^{\ast}\big{(}\hat{\Sigma}_{xy}-\lambda_{n}z^{\lambda_{n}}\big{)}-W^{\ast}\big{(}\hat{\Sigma}_{xy}\big{)}\Big{)}\leqslant\frac{1}{\lambda_{n}}\Big{(}W^{\ast}\big{(}\hat{\Sigma}_{xy}-\lambda_{n}z^{0}\big{)}-W^{\ast}\big{(}\hat{\Sigma}_{xy}\big{)}\Big{)}$
$\displaystyle=$ $\displaystyle\langle-z^{0},\nabla
W^{\ast}(\hat{\Sigma}_{xy})\rangle+\frac{\lambda_{n}}{2}\big{\langle}z^{0},\nabla^{2}W^{\ast}(\hat{\Sigma}_{xy})z^{0}\big{\rangle}+R(\hat{\Sigma}_{xy}+\lambda_{n}z^{0})\lambda^{2}_{n}$
$\displaystyle\leqslant$ $\displaystyle\langle-z^{\lambda_{n}},\nabla
W^{\ast}(\hat{\Sigma}_{xy})\rangle+\frac{\lambda_{n}}{2}\big{\langle}z^{0},\nabla^{2}W^{\ast}(\hat{\Sigma}_{xy})z^{0}\big{\rangle}+R(\hat{\Sigma}_{xy}+\lambda_{n}z^{0})\lambda^{2}_{n},$

where we used the optimality of $z^{0}$ and of $z^{\lambda_{n}}$. We conclude
that

$\displaystyle\frac{1}{2}\big{\langle}z^{\lambda_{n}},\nabla^{2}W^{\ast}(\hat{\Sigma}_{xy})z^{\lambda_{n}}\big{\rangle}+R(\hat{\Sigma}_{xy}+\lambda_{n}z^{\lambda_{n}})\lambda_{n}\leqslant\frac{1}{2}\big{\langle}z^{0},\nabla^{2}W^{\ast}(\hat{\Sigma}_{xy})z^{0})\big{\rangle}+R(\hat{\Sigma}_{xy}+\lambda_{n}z^{0})\lambda_{n}.$

Taking the limit establishes 17. Since $z^{0}$ was an arbitrary element in
$\partial\left\|{\widehat{A}}\right\|_{1}$, we obtain that the limit of
$z^{\lambda_{n}}$ is

$\displaystyle
z^{\ast}=\operatorname*{arg\,min}_{z}\Big{\langle}z,\Big{(}\nabla^{2}W({\widehat{A}})\Big{)}^{-1}z\Big{\rangle}\quad\text{subject
to}\quad z\in\partial\left\|{\widehat{A}}\right\|_{1}$

where we used 15. Finally, observe that $z^{\ast}$ was an arbitrary limit
point of $z^{\lambda}$ and we showed that all limit points are the same; this
is enough to conclude the result. ∎

##  Appendix B Proof of Proposition 7

The proof of this statement relies on strong convexity of $\mathcal{J}_{n}$.
Similar results have been proven in the context of entropic optimal transport
(e.g. [Genevay et al., 2019]). Our proof is similar to the approach taken in
[Rigollet and Stromme, 2022].

Let $(A_{\infty},f_{\infty},g_{\infty})$ minimize ($\mathcal{K}_{\infty}$).
Note that $p_{\infty}\alpha\otimes\beta$ with

$p_{\infty}(x,y)=\exp\left(\frac{2}{\varepsilon}\left(\Phi
A_{\infty}(x,y)+f_{\infty}(x)+g_{\infty}(y)\right)\right)$

minimizes ($\mathrm{iOT-}\ell_{1}(\hat{\pi})$). Let

$P_{\infty}=\frac{1}{n^{2}}(p_{\infty}(x_{i},y_{j}))_{i,j},\quad
F_{\infty}=(f_{\infty}(x_{i}))_{i}\quad\text{and}\quad
G_{\infty}=(g_{\infty}(y_{j}))_{j}.$

Note that by optimality of $(A_{\infty},f_{\infty},g_{\infty})$ and
$\pi_{\infty}=p_{\infty}\alpha\otimes\beta$,

$\lambda\left\|A_{\infty}\right\|_{1}+\sup_{\pi\in\mathcal{U}(\alpha,\beta)}\langle\Phi
A_{\infty},\,\pi-\hat{\pi}\rangle-\frac{\varepsilon}{2}\mathrm{KL}(\pi|\alpha\otimes\beta)\leqslant\frac{\varepsilon}{2}$

where we plug in $(A,F,G)=(0,0,0)$ into ($\mathcal{K}_{\infty}$) for the RHS.
Plugging in $\hat{\pi}$ into the LHS, we have

$\lambda\left\|A_{\infty}\right\|_{1}\leqslant\frac{\varepsilon}{2}+\frac{\varepsilon}{2}\mathrm{KL}(\hat{\pi}|\alpha\otimes\beta)\leqslant
C\varepsilon,$

with the implicit constant $C$ being dependent on $\hat{\pi}$. So,
$\left\|A_{\infty}\right\|_{1}=\mathcal{O}(\varepsilon/\lambda)$ and due to
the uniform bounds on $f_{\infty},g_{\infty}$ from Lemma 13, $p_{\infty}$ is
uniformly bounded away from 0 by $\exp(-C/\lambda)$ for some constant $C$ that
depends on $\hat{\pi}$.

Let $P_{n}$ minimise ($\mathcal{P}_{n}$), we know it is of the form

$P_{n}=\frac{1}{n^{2}}\exp\left(\frac{2}{\varepsilon}\left(\Phi_{n}A_{n}+F_{n}\oplus
G_{n}\right)\right).$

for vectors $A_{n},F_{n},G_{n}$.

The ‘certificates’ are $\Phi^{*}\varphi_{\infty}$ and
$z_{n}\triangleq\Phi_{n}^{*}\varphi_{n}$ where

$\varphi_{\infty}=\frac{1}{\lambda}\left(p_{\infty}\alpha\otimes\beta-\hat{\pi}\right)\quad\text{and}\quad\varphi_{n}=\frac{1}{\lambda}\left(P_{n}-\hat{P}_{n}\right).$

Note that
$z_{n}=\Phi^{*}\varphi_{\infty}=-\frac{1}{\lambda}\nabla_{A}\mathcal{L}(A^{\infty},\hat{\pi})$.
The goal is to bound
$\left\|\Phi^{*}\varphi_{\infty}-\Phi_{n}^{*}\varphi_{n}\right\|_{\infty}$ so
that nondegeneracy of $z_{\infty}$ would imply nondegeneracy of
$\Phi_{n}^{*}\varphi_{n}$. Note that
$\Phi_{n}^{*}P_{n}=\Phi^{*}\hat{\pi}_{n}$. So, by the triangle inequality,

$\displaystyle\left\|\Phi^{*}\varphi_{\infty}-\Phi_{n}^{*}\varphi_{n}\right\|_{\infty}$
$\displaystyle\leqslant\frac{1}{\lambda}\left\|\Phi^{*}\left(p_{\infty}\alpha\otimes\beta\right)-\Phi_{n}^{*}P_{n}\right\|_{\infty}+\frac{1}{\lambda}\left\|\Phi^{*}\left(\hat{\pi}_{n}-\hat{\pi}\right)\right\|_{\infty}$
$\displaystyle\leqslant\frac{1}{\lambda}\left\|\Phi_{n}^{*}P_{\infty}-\Phi_{n}^{*}P_{n}\right\|_{\infty}+\frac{1}{\lambda}\left\|\Phi^{*}\left(p_{\infty}\alpha\otimes\beta\right)-\Phi_{n}^{*}P_{\infty}\right\|_{\infty}+\frac{1}{\lambda}\left\|\Phi^{*}\left(\hat{\pi}_{n}-\hat{\pi}\right)\right\|_{\infty}$

The last two terms on the RHS can be controlled using Proposition 21, and are
bounded by $\mathcal{O}(mn^{-1/2})$ with probability at least
$1-\mathcal{O}(\exp(-m^{2}))$. For the first term on the RHS, letting
$Z=P_{\infty}-P_{n}$,

$\displaystyle\left\|\Phi_{n}^{*}Z\right\|_{\infty}$
$\displaystyle=\left\|\sum_{i,j=1}^{n}\mathbf{C}(x_{i},y_{j})Z_{i,j}\right\|_{\infty}$
$\displaystyle\leqslant\left\|\mathbf{C}\right\|_{\infty}\sqrt{\frac{1}{n^{2}}\sum_{i,j}\left(\exp\left(\frac{2}{\varepsilon}(\Phi_{n}A_{n}+F_{n}\oplus
G_{n})\right)-\exp\left(\frac{2}{\varepsilon}(\Phi_{n}A_{\infty}+F_{\infty}\oplus
G_{\infty})\right)\right)_{i,j}^{2}}$

Let $L\triangleq 2(\left\|\Phi_{n}A_{n}+F_{n}\oplus
G_{n})\right\|_{\infty}\vee\left\|\Phi_{n}A_{\infty}+F_{\infty}\oplus
G_{\infty})\right\|_{\infty})$. By Lipschitz continuity of the exponential

$\displaystyle\frac{1}{n^{2}}\sum_{i,j}\left(\exp\left(\frac{2}{\varepsilon}(\Phi_{n}A_{n}+F_{n}\oplus
G_{n})\right)-\exp\left(\frac{2}{\varepsilon}(\Phi_{n}A_{\infty}+F_{\infty}\oplus
G_{\infty})\right)\right)_{i,j}^{2}$ (18)
$\displaystyle\leqslant\frac{4\exp(L/\varepsilon)}{\varepsilon^{2}n^{2}}\sum_{i,j}\left((\Phi_{n}A_{n}+F_{n}\oplus
G_{n})-(\Phi_{n}A_{\infty}+F_{\infty}\oplus G_{\infty})\right)_{i,j}^{2}$ (19)
$\displaystyle\leqslant\frac{12\exp(L/\varepsilon)}{\varepsilon^{2}}\left(\frac{1}{n^{2}}\sum_{i,j}(\Phi_{n}A_{n}-\Phi_{n}A_{\infty})_{i,j}^{2}+\frac{1}{n}\sum_{i}(F_{n}-F_{\infty})_{i}^{2}+\frac{1}{n}\sum_{j}(G_{n}-G_{\infty})_{j}^{2}\right)$
(20)

By strong convexity properties of $\mathcal{J}_{n}$ and hence
$\mathcal{K}_{n}$, it can be shown (see Prop 15) that (20) is upper bounded up
to a constant by

$\displaystyle\varepsilon^{-1}\exp(L/\varepsilon)\left(\mathcal{K}_{n}(F_{\infty},G_{\infty},A_{\infty})-\mathcal{K}_{n}(F_{n},G_{n},A_{n})\right)$
$\displaystyle\leqslant$
$\displaystyle\frac{\exp(L/\varepsilon)}{4}\Bigg{(}\left\|n^{-2}(\Phi^{*}_{n}\Phi_{n})^{-1}\right\|\left\|(\partial_{A}\mathcal{J}_{n}(A_{\infty},F_{\infty},G_{\infty})+\lambda\xi_{\infty})\right\|^{2}$
$\displaystyle\qquad\qquad+n\left\|\partial_{F}\mathcal{J}_{n}(A_{\infty},F_{\infty},G_{\infty})\right\|^{2}+n\left\|\partial_{G}\mathcal{J}_{n}(A_{\infty},F_{\infty},G_{\infty})\right\|^{2}\Bigg{)}$

where
$\xi_{\infty}=\frac{1}{\lambda}(\Phi^{*}\hat{\pi}-\Phi^{*}(p_{\infty}\alpha\otimes\beta))\in\partial\left\|A_{\infty}\right\|_{1}$.
By Lemma 12 and Lemma 13, $L=\mathcal{O}(\varepsilon/\lambda)$ with
probability at least $1-\mathcal{O}(\exp(-m^{2}))$ if $\lambda\gtrsim
mn^{-\frac{1}{2}}$. Finally,

$\displaystyle\left\|(\partial_{A}\mathcal{J}_{n}(A_{\infty},F_{\infty},G_{\infty})+\lambda\xi_{\infty})\right\|^{2}=\left\|\Phi_{n}^{*}\hat{P}+\Phi_{n}^{*}P_{\infty}+\lambda\xi_{\infty}\right\|^{2}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\leqslant
2\left(\left\|\Phi_{n}^{*}\hat{P}-\Phi^{*}\hat{\pi}\right\|^{2}+\left\|\Phi_{n}^{*}P_{\infty}-\Phi^{*}(p_{\infty}\alpha\otimes\beta)\right\|^{2}\right)$
$\displaystyle
n\left\|\partial_{F}\mathcal{J}_{n}(A_{\infty},F_{\infty},G_{\infty})\right\|^{2}_{2}=\frac{1}{n}\sum_{i=1}^{n}\left(1-\frac{1}{n}\sum_{j=1}^{n}p_{\infty}(x_{i},y_{j})\right)^{2}$
$\displaystyle
n\left\|\partial_{G}\mathcal{J}_{n}(A_{\infty},F_{\infty},G_{\infty})\right\|^{2}_{2}\Bigg{)}=\frac{1}{n}\sum_{j=1}^{n}\left(1-\frac{1}{n}\sum_{i=1}^{n}p_{\infty}(x_{i},y_{j})\right)^{2}\Bigg{)}.$

We show in Propositions 17, 20 and 21 that these are bounded by
$\mathcal{O}(m^{2}n^{-1})$ with probability at least
$1-\mathcal{O}(\exp(-m^{2}))$ and from Proposition 24, assuming
$n\gtrsim\log(2s)$, $\left\|(n^{-2}\Phi_{n}^{*}\Phi_{n})^{-1}\right\|\lesssim
1$ with probability at least $1-\mathcal{O}(\exp(-n))$.

So,
$\lambda\left\|\Phi^{*}\varphi_{\infty}-\Phi_{n}^{*}\varphi_{n}\right\|_{\infty}\lesssim\exp(C/\lambda)\frac{m}{\sqrt{n}}$
with probability at least $1-\exp(-m^{2})$. The second statement follows by
combining our bound for (20) with the fact that
$\left\|(n^{-2}\Phi_{n}^{*}\Phi_{n})^{-1}\right\|\lesssim 1$.

In the following subsections, we complete the proof by establishing the
required strong convexity properties of $\mathcal{J}_{n}$ and bound
$\nabla\mathcal{J}_{n}$ at $(A_{\infty},F_{\infty},G_{\infty})$ using
concentration inequalities. The proofs to some of these results are verbatim
to the results of [Rigollet and Stromme, 2022] for deriving sampling
complexity bounds in eOT, although they are included due to the difference in
our setup.

###  B.1 Strong convexity properties of $\mathcal{J}_{n}$

In this section, we present some of the key properties of $\mathcal{J}_{n}$.

Recall that since we assume that $\alpha,\beta$ are compactly supported, up to
a rescaling of the space, we assume without loss of generality that for all
$k$, $\left\lvert\mathbf{C}_{k}(x,y)\right\rvert\leqslant 1$.

######  Lemma 12.

Let $A,F,G$ minimize $\mathcal{K}_{n}$. Let $\lambda\gtrsim mn^{-\frac{1}{2}}$
for some $m>0$. Then, with probability at least $1-\mathcal{O}(\exp(-m^{2}))$

$\left\|A\right\|_{1}\lesssim(C+1)\varepsilon/\lambda$

where the constant $C>0$ depends on $d\hat{\pi}/d(\alpha\otimes\beta)$.

###### Proof.

Recall $\hat{P}_{n}=\frac{1}{n}\mathrm{Id}$. By optimality of $A,F,G$,

$\lambda\left\|A\right\|_{1}-\langle\hat{P}_{n},\,\Phi_{n}A\rangle+\langle
P,\,\Phi_{n}A\rangle-\frac{\varepsilon}{2}H(P)\leqslant\frac{\varepsilon}{2}.$
(21)

for all $P$ satisfying the marginal constraints
$P\mathds{1}=\frac{1}{n}\mathds{1}$ and
$P^{\top}\mathds{1}=\frac{1}{n}\mathds{1}$. By plugging in $P=\hat{P}_{n}$, we
have $H(P)=\frac{1}{n}\sum_{i=1}^{n}\log(n)-1=\log(n)-1$. So,

$\lambda\left\|A\right\|_{1}\leqslant\varepsilon\log(n).$

To obtain a bound that does not grow with $n$, the idea is to take
$Q_{n}=\frac{1}{n^{2}}(\hat{p}(x_{i},y_{j}))_{i,j\in[n]}$ where
$\hat{p}=d\hat{\pi}/d(\alpha\otimes\beta)$. Note that $Q_{n}$ will only
approximately satisfy the marginal constraints
$Q_{n}\mathds{1}\approx\frac{1}{n}\mathds{1}$ and
$Q_{n}^{\top}\mathds{1}\approx\frac{1}{n}\mathds{1}$ (this approximation can
be made precise using Proposition 17). So, we insert into the above inequality
$P=\tilde{Q}_{n}$ with $\tilde{Q}_{n}$ being the projection of $Q_{n}$ onto
the constraint set
$\mathcal{U}(\frac{1}{n}\mathds{1}_{n},\frac{1}{n}\mathds{1}_{n})$.

By [Altschuler et al., 2017, Lemma 7], the projection $\tilde{Q}_{n}$
satisfies

$\left\|\tilde{Q}_{n}-Q_{n}\right\|_{1}\leqslant
2\left\|Q_{n}\mathds{1}-\frac{1}{n}\mathds{1}\right\|_{1}+2\left\|Q_{n}^{\top}\mathds{1}-\frac{1}{n}\mathds{1}\right\|_{1}.$

By Proposition 17, with probability at least $1-\mathcal{O}(\exp(-m^{2}))$,

$\left\|Q_{n}\mathds{1}-\frac{1}{n}\mathds{1}\right\|_{1}=\frac{1}{n}\sum_{i=1}^{n}\left\lvert\frac{1}{n}\sum_{j=1}^{n}\hat{p}(x_{i},y_{j})-1\right\rvert\leqslant\sqrt{\frac{1}{n}\sum_{i=1}^{n}\left\lvert\frac{1}{n}\sum_{j=1}^{n}\hat{p}(x_{i},y_{j})-1\right\rvert^{2}}=\mathcal{O}\left(\frac{m}{\sqrt{n}}\right).$

Similarly,
$\left\|Q_{n}\mathds{1}-\frac{1}{n}\mathds{1}\right\|_{1}=\mathcal{O}(\log(n)n^{-\frac{1}{2}})$
with high probability. Moreover, since $\hat{p}$ is uniformly bounded from
below, one can check from the projection procedure of [Altschuler et al.,
2017] that $\tilde{Q}_{n}$ is also uniformly bounded from above and away from
zero with constants that dependent only on $\hat{p}$. So,
$H(\tilde{Q}_{n})=C$, where $C$ is a constant depends on the lower bound of
$d\hat{\pi}/d(\alpha\otimes\beta)$.

Plugging in $\tilde{Q}_{n}$ into (21), we have

$\lambda\left\|A\right\|_{1}-\langle\hat{P}_{n}-\tilde{Q}_{n},\,\Phi_{n}A\rangle\lesssim
C\varepsilon+\varepsilon.$

Note that

$\displaystyle\left\lvert\langle\tilde{Q}_{n}-Q_{n},\,\Phi_{n}A\rangle\right\rvert=\left\lvert\sum_{k=1}^{s}\sum_{i,j=1}^{n}A_{k}(C_{k})_{i,j}(\tilde{Q}_{n}-Q_{n})_{i,j}\right\rvert$
(22)
$\displaystyle\leqslant\max_{k}\left\|C_{k}\right\|_{\infty}\left\|A\right\|_{1}\left\|Q_{n}-\tilde{Q}_{n}\right\|_{1}$
(23)

Moreover,

$\displaystyle\langle\hat{P}_{n}-Q_{n},\,\Phi_{n}A\rangle$
$\displaystyle\leqslant\left\|A\right\|\left\|\Phi_{n}^{*}(\hat{P}_{n}-Q_{n})\right\|$
(24)

By Proposition 20 and 21, with probability at least
$1-\mathcal{O}(\exp(-m^{2}))$,

$\mathbb{E}\left\|\Phi_{n}^{*}\hat{P}_{n}-\Phi^{*}\hat{\pi}\right\|=\mathcal{O}(mn^{-\frac{1}{2}})\quad\text{and}\quad\mathbb{E}\left\|\Phi_{n}^{*}Q_{n}-\Phi^{*}\hat{\pi}\right\|=\mathcal{O}(mn^{-\frac{1}{2}}).$

So, we have

$\left(\lambda-\left\|\Phi_{n}^{*}(\hat{P}_{n}-Q_{n})\right\|-\left\|\tilde{Q}_{n}-Q_{n}\right\|_{1}\right)\left\|A\right\|_{1}\lesssim(C+1)\varepsilon,$

and hence, with probability at least $1-\mathcal{O}(\exp(-m^{2}))$

$\left\|A\right\|_{1}\lesssim(C+1)\varepsilon/\lambda$

if $\lambda\gtrsim mn^{-\frac{1}{2}}$.

∎

######  Lemma 13.

Let $(A,F,G)$ minimize ($\mathcal{K}_{n}$). Let $C_{A}=\Phi_{n}A$. Then,

$F_{i}\in[-3\left\|C_{A}\right\|_{\infty},\left\|C_{A}\right\|_{\infty}]\quad\text{and}\quad
G_{i}\in[-2\left\|C_{A}\right\|_{\infty},2\left\|C_{A}\right\|_{\infty}].$

Moreover,
$\exp(4\left\|C_{A}\right\|_{\infty})\geqslant\exp\left(\frac{2}{\varepsilon}\left(F_{i}+G_{j}+(\Phi_{n}A)_{i,j}\right)\right)\geqslant\exp(-6\left\|C_{A}\right\|_{\infty})$.
Note that $\left\|C_{A}\right\|_{\infty}\leqslant\left\|A\right\|_{1}$.

If $(A,f,g)$ minimize ($\mathcal{K}_{\infty}$). Let $c_{A}=\Phi A$. Then, for
all $x,y$

$f(x)\in[-3\left\|c_{A}\right\|_{\infty},\left\|c_{A}\right\|_{\infty}]\quad\text{and}\quad
g(y)\in[-2\left\|c_{A}\right\|_{\infty},2\left\|c_{A}\right\|_{\infty}].$

Moreover,
$\exp(4\left\|c_{A}\right\|_{\infty})\geqslant\exp\left(\frac{2}{\varepsilon}\left(f\oplus+g+(\Phi
A)_{i,j}\right)\right)\geqslant\exp(-6\left\|c_{A}\right\|_{\infty})$.

###### Proof.

This proof is nearly identical to [Rigollet and Stromme, 2022, Prop. 10]: Let
$A,F,G$ minimize ($\mathcal{K}_{n}$).

By the marginal constraints for
$P_{n}\triangleq\frac{1}{n^{2}}\left(\exp\left(\frac{2}{\varepsilon}(F\oplus
G+\Phi_{n}A)\right)\right)_{i,j}$ given in ($\mathcal{P}_{n}$),

$\begin{split}1&=\frac{1}{n}\sum_{j}\exp\left(\frac{2}{\varepsilon}(\Phi
A(x_{i},y_{j})+F_{i}+G_{j})\right)\\\
&\geqslant\exp\left(\frac{2}{\varepsilon}(-\left\|C_{A}\right\|_{\infty}+F_{i})\right)\sum_{j}\frac{1}{n}\exp\left(2G_{j}/\varepsilon\right)\geqslant\exp\left(\frac{2}{\varepsilon}(-\left\|C_{A}\right\|_{\infty}+F_{i})\right)\end{split}$
(25)

where we use Jensen’s inequality and the assumption that $\sum_{j}G_{j}=0$ for
the second. So, $F_{i}\leqslant\left\|C_{A}\right\|_{\infty}$ for all
$i\in[n]$. Using the other marginal constraint for $P_{n}$ along with this
bound on $F_{i}$ implies that

$1=\frac{1}{n}\sum_{i}\exp\left(\frac{2}{\varepsilon}(\Phi
A(x_{i},y_{j})+F_{i}+G_{j})\right)\leqslant\exp\left(\frac{2}{\varepsilon}\left(2\left\|C_{A}\right\|_{\infty}+G_{j})\right)\right)$

So, $G_{n,j}\geqslant-2\left\|C_{A}\right\|_{\infty}$.

To prove the reverse bounds, we now show that $\sum_{i}F_{i}$ is lower
bounded: Note that since $H(P)\geqslant-1$, by duality between
($\mathcal{P}_{n}$) and ($\mathcal{K}_{n}$),

$\left\|C_{A}\right\|_{\infty}+\frac{\varepsilon}{2}\geqslant\langle\Phi_{n}A,\,P\rangle-\frac{\varepsilon}{2}H(P)=-\frac{1}{n}\sum_{i}F_{i}-\frac{1}{n}\sum_{j}G_{j}+\frac{\varepsilon}{2n^{2}}\sum_{i,j}\exp\left(\frac{2}{\varepsilon}(F_{i}+G_{j}+(\Phi_{n}A)_{ij}\right)$

By assumption, $\sum_{j}G_{j}=0$ and
$\sum_{i,j}\exp\left(\frac{2}{\varepsilon}(F_{i}+G_{j}+(\Phi_{n}A)_{ij}\right)=n^{2}$.
So,

$\frac{1}{n}\sum_{i}F_{i}\geqslant-\left\|C_{A}\right\|_{\infty}.$

By repeating the argument in (25), we see that

$1\geqslant\exp(2/\varepsilon(-\left\|C_{A}\right\|_{\infty}+G_{j}))\exp\left(\frac{2}{n\varepsilon}\sum_{i}F_{i}\right)\geqslant\exp(2/\varepsilon(-2\left\|C_{A}\right\|_{\infty}+G_{j})).$

So, $G_{j}\leqslant 2\left\|C_{A}\right\|_{\infty}$ and
$F_{i}\geqslant-3\left\|C_{A}\right\|_{\infty}$.

The proof for ($\mathcal{K}_{\infty}$) is nearly identical and hence omitted.
∎

Similarly to [Rigollet and Stromme, 2022, Lemma 11], we derive the following
strong convexity bound for $\mathcal{J}_{n}$:

######  Lemma 14.

The functional $\mathcal{J}_{n}$ is strongly convex with

$\begin{split}\mathcal{J}_{n}(A^{\prime},F^{\prime},G^{\prime})\geqslant\mathcal{J}_{n}(A,F,G)&+\langle\nabla\mathcal{J}_{n}(A,F,G),\,(A^{\prime},F^{\prime},G^{\prime})-(A,F,G)\rangle\\\
&+\frac{\exp(-L/\varepsilon)}{\varepsilon}\left(\frac{1}{n^{2}}\left\|\Phi_{n}(A-A^{\prime})\right\|^{2}_{2}+\frac{1}{n}\left\|F-F^{\prime}\right\|^{2}_{2}+\frac{1}{n}\left\|G-G^{\prime}\right\|_{2}^{2}\right),\end{split}$
(26)

for some
$L=\mathcal{O}(\left\|A\right\|_{1}\vee\left\|A^{\prime}\right\|_{1})$.

###### Proof.

To establish the strong convexity inequality, let

$h(t)\triangleq\mathcal{J}((1-t)A+tA^{\prime},(1-t)f+tf^{\prime},(1-t)g+tg^{\prime}).$

It suffices to find $\delta>0$ such that for all $t\in[0,1]$,

$h^{\prime\prime}(t)\geqslant\delta\left(\frac{1}{n^{2}}\left\|\Phi_{n}(A-A^{\prime})\right\|^{2}_{2}+\frac{1}{n}\left\|F-F^{\prime}\right\|^{2}_{2}+\frac{1}{n}\left\|G-G^{\prime}\right\|_{2}^{2}\right).$
(27)

Let $Z_{t}\triangleq((1-t)F+tF^{\prime})\oplus((1-t)G+tG^{\prime})+((1-t)\Phi
A+t\Phi A^{\prime})$. Note that

$\displaystyle h^{\prime\prime}(t)=\frac{2}{\varepsilon
n^{2}}\sum_{i,j}\exp\left(\frac{2}{\varepsilon}(Z_{t})_{i,j}\right)(F^{\prime}_{i}-F_{i}+G^{\prime}_{i}-G_{i}+(\Phi
A^{\prime})_{i,j}-(\Phi A)_{i,j})^{2}$ (28)

Since
$\left\|c_{A}\right\|_{\infty}\vee\left\|c_{A^{\prime}}\right\|_{\infty}\leqslant
L$, by Lemma 13,
$\left\|Z_{t}\right\|_{\infty}\lesssim\left\|c_{A}\right\|_{\infty}\vee\left\|c_{A^{\prime}}\right\|_{\infty}$.
So,

$h^{\prime\prime}(t)\geqslant\frac{2}{\varepsilon
n^{2}}\exp(-L/\varepsilon)\sum_{i,j}(F^{\prime}_{i}-F_{i}+G^{\prime}_{i}-G_{i}+(\Phi
A^{\prime})_{i,j}-(\Phi A)_{i,j})^{2}.$

By expanding out the brackets and using $\sum_{i}G_{i}=0$ and since $C_{k}$
are centred ( $\sum_{i}(C_{k})_{i,j}=0$ and $\sum_{j}(C_{k})_{i,j}=0$), (27)
holds with $\delta=\frac{2}{\varepsilon}\exp(-L/\varepsilon)$. ∎

Based on this strong convexity result, we have the following bound.

######  Proposition 15.

Let $(A_{n},F_{n},G_{n})$ minimise $\mathcal{F}_{n}$, then for all
$(A,F,G)\in\mathcal{S}$ such that $n^{-2}\exp(F\oplus G+\Phi_{n}A)$ satisfy
the marginal constraints of ($\mathcal{P}_{n}$),

$\mathcal{K}_{n}(A,F,G)-\mathcal{K}_{n}(A_{n},F_{n},G_{n})\geqslant\frac{\exp(-L/\varepsilon)}{\varepsilon}\left(\frac{1}{n^{2}}\left\|\Phi_{n}(A-A_{n})\right\|^{2}_{2}+\frac{1}{n}\left\|F-F_{n}\right\|^{2}_{2}+\frac{1}{n}\left\|G-G_{n}^{\prime}\right\|_{2}^{2}\right).$

and

$\displaystyle\mathcal{K}_{n}(A,F,G)-\mathcal{K}_{n}(A_{n},F_{n},G_{n})\leqslant$
$\displaystyle\frac{\varepsilon\exp(L/\varepsilon)}{4}\Bigg{(}\left\|n^{-2}(\Phi^{*}_{n}\Phi_{n})^{-1}\right\|\left\|(\partial_{A}\mathcal{J}_{n}(A,F,G)+\lambda\xi)\right\|^{2}$
$\displaystyle+n\left\|\partial_{F}\mathcal{J}_{n}(A,F,G)\right\|^{2}_{2}+n\left\|\partial_{G}\mathcal{J}_{n}(A,F,G)\right\|^{2}_{2}\Bigg{)}.$

where $L=\mathcal{O}(\left\|A\right\|_{1}\vee\left\|A_{n}\right\|_{1})$.

###### Proof.

By strong convexity of $\mathcal{J}_{n}$, we can show that for any
$(A,f,g),(A^{\prime},f^{\prime},g^{\prime})\in\mathcal{S}$ and any
$\xi\in\partial\left\|A\right\|_{1}$,

$\begin{split}\mathcal{K}_{n}(A^{\prime},F^{\prime},G^{\prime})\geqslant\mathcal{K}_{n}(A,F,G)&+\langle\nabla\mathcal{J}_{n}(A,F,G),\,(A^{\prime},F^{\prime},G^{\prime})-(A,F,G)\rangle+\lambda\langle\xi,\,A^{\prime}-A\rangle\\\
&+\frac{\delta}{2}\left(\frac{1}{n^{2}}\left\|\Phi_{n}(A-A^{\prime})\right\|^{2}_{2}+\frac{1}{n}\left\|F-F^{\prime}\right\|^{2}_{2}+\frac{1}{n}\left\|G-G^{\prime}\right\|_{2}^{2}\right),\end{split}$
(29)

where $\delta=\frac{2\exp(-L/\varepsilon)}{\varepsilon}$ with
$L=\mathcal{O}(\left\|A\right\|_{1}\vee\left\|A^{\prime}\right\|_{1})$. The
first statement follows by letting $(A,F,G)=(A_{n},F_{n},G_{n})$ in the above
inequality.

To prove the second statement, let

$\displaystyle M\triangleq$
$\displaystyle\langle\nabla\mathcal{J}_{n}(A,F,G),\,(A^{\prime},F^{\prime},G^{\prime})-(A,F,G)\rangle+\lambda\langle\xi,\,A^{\prime}-A\rangle$
$\displaystyle+\frac{\delta}{2}\left(\frac{1}{n^{2}}\left\|\Phi_{n}(A-A^{\prime})\right\|^{2}_{2}+\frac{1}{n}\left\|F-F^{\prime}\right\|^{2}_{2}+\frac{1}{n}\left\|G-G^{\prime}\right\|_{2}^{2}\right).$

By minimising over $(A^{\prime},F^{\prime},G^{\prime})$, note that

$M\geqslant-\frac{1}{2\delta}\left(n^{2}\left\|(\partial_{A}\mathcal{J}_{n}(A,F,G)+\lambda\xi)\right\|^{2}_{(\Phi_{n}^{*}\Phi_{n})^{-1}}+n\left\|\partial_{F}\mathcal{J}_{n}(A,F,G)\right\|^{2}_{2}+n\left\|\partial_{G}\mathcal{J}_{n}(A,F,G)\right\|^{2}_{2}\right),$

So,

$-M\leqslant\frac{1}{2\delta}\left(\left\|n^{-2}(\Phi^{*}_{n}\Phi_{n})^{-1}\right\|\left\|(\partial_{A}\mathcal{J}_{n}(A,F,G)+\lambda\xi)\right\|^{2}+n\left\|\partial_{F}\mathcal{J}_{n}(A,F,G)\right\|^{2}_{2}+n\left\|\partial_{G}\mathcal{J}_{n}(A,F,G)\right\|^{2}_{2}\right)$

Finally, note that
$\mathcal{K}_{n}(A,F,G)-\mathcal{K}_{n}(A^{\prime},F^{\prime},G^{\prime})\geqslant\mathcal{K}_{n}(A,F,G)-\mathcal{K}_{n}(A_{n},F_{n},G_{n})$
by optimality of $A_{n},F_{n},G_{n}$. ∎

###  B.2 Concentration bounds

######  Lemma 16 (McDiarmid’s inequality).

Let $f:\mathcal{X}^{n}\to\mathbb{R}$ satisfy the bounded differences property:

$\sup_{x_{i}^{\prime}\in\mathcal{X}}\left\lvert
f(x_{1},\ldots,x_{i-1},x_{i},x_{i+1},\ldots,x_{n})-f(x_{1},\ldots,x_{i-1},x_{i}^{\prime},x_{i+1},\ldots,x_{n})\right\rvert\leqslant
c.$

Then, given $X_{1},\ldots,X_{n}$ random variables with $X_{i}\in\mathcal{X}$,
for any $t>0$,

$\mathbb{P}(.\left\lvert
f(X_{1},\ldots,X_{n})-\mathbb{E}[f(X_{1},\ldots,X_{n})]\right\rvert\geqslant
t)\leqslant 2\exp(-2t^{2}/(nc^{2})).$

Given random vectors $X$ and $Y$, denote $\mathrm{Cov}(X,Y)=\mathbb{E}\langle
Y-\mathbb{E}[Y],\,X-\mathbb{E}[X]\rangle$ and
$\mathrm{Var}(X)=\mathrm{Cov}(X,X)$.

######  Proposition 17.

Let $\pi$ have marginals $\alpha$ and $\beta$ and let
$p=\frac{d\pi}{d(\alpha\otimes\beta)}$. Let $(x_{i},y_{i})\sim\hat{\pi}$ where
$\hat{\pi}$ has marginals $\alpha,\beta$. Assume that
$\left\|p\right\|_{\infty}\leqslant b$. Then,

$\displaystyle\mathbb{E}\left[\frac{1}{n}\sum_{j=1}^{n}\left(1-\frac{1}{n}\sum_{i=1}^{n}p(x_{i},y_{j})\right)^{2}\right]\leqslant\frac{(b+1)^{2}}{n}$

and

$\frac{1}{n}\sum_{j=1}^{n}\left(1-\frac{1}{n}\sum_{i=1}^{n}p(x_{i},y_{j})\right)^{2}\leqslant\frac{(t+b+1)^{2}}{n}.$

with probability at least $1-\exp(-t^{2}/(4b^{2}))$.

######  Remark 18.

The bounds also translate to an $\ell_{1}$ norm on the marginal errors, since
by Cauchy-Schwarz,

$\displaystyle\frac{1}{n}\sum_{j=1}^{n}\left\lvert
1-\frac{1}{n}\sum_{i=1}^{n}p(x_{i},y_{j})\right\rvert$
$\displaystyle\leqslant\sqrt{\sum_{j=1}^{n}\frac{1}{n}\left(1-\frac{1}{n}\sum_{i=1}^{n}p(x_{i},y_{j})\right)^{2}}.$

Moreover, by Jensen’s inequality
$\mathbb{E}\frac{1}{n}\sum_{j=1}^{n}\left\lvert
1-\frac{1}{n}\sum_{i=1}^{n}p(x_{i},y_{j})\right\rvert\leqslant\sqrt{\mathbb{E}\sum_{j=1}^{n}\frac{1}{n}\left(1-\frac{1}{n}\sum_{i=1}^{n}p(x_{i},y_{j})\right)^{2}}.$

###### Proof.

$\displaystyle\mathbb{E}\left[\frac{1}{n}\sum_{j=1}^{n}\left(1-\frac{1}{n}\sum_{i=1}^{n}p(x_{i},y_{j})\right)^{2}\right]=\frac{1}{n^{3}}\sum_{ij,k=1}^{n}\mathbb{E}\left((1-p(x_{i},y_{j}))(1-p(x_{k},y_{j}))\right).$

For each $j\in[n]$, we have the following cases for
$u_{j}\triangleq\mathbb{E}\left((1-p(x_{i},y_{j}))(1-p(x_{k},y_{j}))\right)$:

  1. 1.

$i=k=j$, then $u_{j}=\mathbb{E}_{(x,y)\sim\hat{\pi}}(p_{\infty}(x,y)-1)^{2}$.
There is 1 such term.

  2. 2.

$i=j$ and $k\neq j$, then
$u_{j}=\mathbb{E}_{(x,y)\sim\hat{\pi},z\sim\alpha}\left((1-p(x,y))(1-p(z,y))\right)=0$.
There are $n-1$ such terms.

  3. 3.

$i\neq j$ and $k=j$, then
$u_{j}=\mathbb{E}_{(z,y)\sim\hat{\pi},x\sim\alpha}\left((1-p(x,y))(1-p(z,y))\right)=0$.
There are $n-1$ such terms.

  4. 4.

$i=k$ and $i\neq j$, then
$u_{j}=\mathbb{E}_{x\sim\alpha,y\sim\beta}(1-p(x,y))^{2}$ and there are
$(n-1)$ such terms.

  5. 5.

$i,j,k$ all distinct. Then, $u_{j}=0$ and there are $(n-1)(n-2)$ such terms.

Therefore,

$\frac{1}{n^{3}}\sum_{ij,k=1}^{n}\mathbb{E}\left((1-p(x_{i},y_{j}))(1-p(x_{k},y_{j}))\right)=\frac{1}{n^{2}}\mathbb{E}_{(x,y)\sim\hat{\pi}}(p(x,y)-1)^{2}+\frac{n-1}{n^{2}}\mathbb{E}_{x\sim\alpha,y\sim\beta}(1-p(x,y))^{2}$

Using $\left\lvert 1-p(x,y)\right\rvert\leqslant b+1$ gives the first
inequality.

Note also that letting
$V=\frac{1}{n}\sum_{i=1}^{n}\left((1-p(x_{i},y_{j}))\right)_{j}$,
$\left\|V\right\|^{2}=\sum_{j=1}^{n}\left(\frac{1}{n}\sum_{i=1}^{n}(1-p(x_{i},y_{j}))\right)^{2}$.
We will apply Lemma 16 to $f(z_{1},\ldots,z_{n})=\left\|V\right\|$. Let
$V^{\prime}=\frac{1}{n}\sum_{i=1}^{n}\left((1-p(x^{\prime}_{i},y^{\prime}_{j}))\right)_{j}$
where $x_{i},y_{j}=x_{j}^{\prime},y_{j}^{\prime}$ for $i,j\geqslant 2$. Then,
for all vectors $u$ of norm 1,

$\displaystyle\langle u,\,V^{\prime}-V\rangle$
$\displaystyle=\sum_{j=1}^{n}u_{j}\frac{1}{n}\sum_{i=1}^{n}(p(x_{i}^{\prime},y_{j}^{\prime})-p(x_{i},y_{j}))$
$\displaystyle=u_{1}\frac{1}{n}\sum_{i=1}^{n}(p(x_{i}^{\prime},y_{1}^{\prime})-p(x_{i},y_{1}))+\sum_{j=2}^{n}u_{j}\frac{1}{n}(p(x_{1}^{\prime},y_{j}^{\prime})-p(x_{1},y_{j}))$
$\displaystyle\leqslant\frac{2b}{n}\sum_{j=1}^{n}u_{j}\leqslant\frac{2b}{\sqrt{n}}.$

So, by the reverse triangle inequality,
$\left\lvert\left\|V\right\|-\left\|V^{\prime}\right\|\right\rvert\leqslant\frac{2b}{\sqrt{n}}$.
It follows that

$n^{-1/2}\left\|V\right\|\leqslant
n^{-1/2}t+n^{-1/2}\mathbb{E}\left\|V\right\|\leqslant
n^{-1/2}t+\sqrt{n^{-1}\mathbb{E}\left\|V\right\|^{2}}\leqslant
n^{-1/2}(t+b+1).$

with probability at least $1-\exp(-t^{2}/(4b^{2}))$ as required.

∎

####  B.2.1 Bounds for the cost parametrization

In the following two propositions (Prop 21 and 20), we assume that $\Phi_{n}$
is defined with $C_{k}=(\mathbf{C}_{k}(x_{i},y_{j}))_{i,j}$ and discuss in the
remark afterwards how to account for the fact that our cost $C_{k}$ in
($\mathcal{P}_{n}$) are centred.

Note that the following classical result is a direct consequence of Lemma 16

######  Lemma 19.

Suppose $Z_{1},\ldots,Z_{m}$ are independent mean zero random variables taking
values in Hilbert space $\mathcal{H}$. Suppose there is some $C>0$ such that
for all $k$, $\left\|Z_{k}\right\|\leqslant C$. Then, for all $t>0$, with
probability at least $1-2\exp(-t)$,

$\left\|\frac{1}{m}\sum_{k=1}^{m}Z_{k}\right\|^{2}\leqslant\frac{8C^{2}t}{m}.$

######  Proposition 20.

Let $C=\max_{x,y}\left\|\mathbf{C}(x,y)\right\|$. Let
$(x_{i},y_{i})_{i=1}^{n}$ be iid drawn from $\hat{\pi}$. Then,

$\left\|\Phi_{n}^{*}\hat{P}-\Phi^{*}\hat{\pi}\right\|\leqslant\frac{\sqrt{8}Ct}{n}$

with probability at least $1-2\exp(-t^{2})$.

###### Proof.

Direct consequence of Lemma 19. ∎

######  Proposition 21.

Let $(x_{i},y_{i})_{i=1}^{n}$ be iid drawn from $\hat{\pi}$, which has
marginals $\alpha,\beta$. Let $P=\frac{1}{n^{2}}(p(x_{i},y_{j}))_{i,j=1}^{n}$
where $\pi$ has marginals $\alpha,\beta$ and
$p=\frac{d\pi}{d(\alpha\otimes\beta)}$. Then,
$\mathbb{E}\left\|\Phi_{n}^{*}P-\Phi^{*}(p\alpha\otimes\beta)\right\|^{2}=\mathcal{O}(n^{-1})$
and

$\left\|\Phi_{n}^{*}P-\Phi^{*}p\alpha\otimes\beta\right\|\lesssim\frac{1+m}{\sqrt{n}}$

with probability at least $1-2\exp(-2m^{2}/(64\left\|p\right\|_{\infty}^{2}))$

###### Proof.

Let $h(x,y)\triangleq p(x,y)\mathbf{C}(x,y)$, then

$\displaystyle\Phi_{n}^{*}P=\frac{1}{n^{2}}\sum_{j=1}^{n}\left(h(x_{j},y_{j})+\sum_{\begin{subarray}{c}i=1\\\
i\neq j\end{subarray}}^{n}h(x_{i},y_{j})\right)$

and

$\mathbb{E}[\Phi_{n}^{*}P-\Phi^{*}p\alpha\otimes\beta]=\frac{1}{n}A^{*}\left(p\hat{\pi}-p\alpha\otimes\beta\right)$

It follows that

$\displaystyle\mathbb{E}\left\|\Phi_{n}^{*}P-\Phi^{*}p\alpha\otimes\beta\right\|^{2}=\frac{1}{n^{2}}\left\|\Phi^{*}\left(p\hat{\pi}-p\alpha\otimes\beta\right)\right\|^{2}+\mathrm{Var}(\Phi_{n}^{*}P)$

and

$\displaystyle\mathrm{Var}(\Phi_{n}^{*}P)=\frac{1}{n^{4}}\sum_{i,j,k,\ell=1}^{n}\mathrm{Cov}(h(x_{i},y_{k}),h(x_{\ell},y_{k}))$

Note that $\mathrm{Cov}(h(x_{i},y_{k}),h(x_{\ell},y_{k}))=0$ whenever
$i,j,k,\ell$ are distinct and there are
$n(n-1)(n-2)(n-3)=n^{4}-6n^{3}+12n^{2}-6n$ such terms, i.e. there are
$6n^{3}-12n^{2}+6n$ nonzero terms in the sum. It follows that
$\mathrm{Var}(\Phi_{n}^{*}P)=\mathcal{O}(n^{-1})$ if
$\left\|\mathbf{C}\right\|$ and $p$ are uniformly bounded.

To conclude, we apply Lemma 16 with
$f(z_{1},\ldots,z_{n})=\left\|\Phi_{n}^{*}P-\Phi^{*}p\alpha\otimes\beta\right\|$
and $z_{i}=(x_{i},y_{i})$. Let
$X(z_{1},\ldots,z_{n})=\Phi_{n}^{*}P-\Phi^{*}p\alpha\otimes\beta$. Then, Note
that for an arbitrary vector $u$,

$\displaystyle\langle
u,\,X(z_{1},z_{2},\ldots,z_{n})-X(z_{1}^{\prime},z_{2},\ldots,z_{n})\rangle=\sum_{i,j}\left(\sum_{k=1}^{s}u_{k}(h(x_{i},y_{j})-h(x_{i}^{\prime},y_{j}^{\prime}))\right)$
$\displaystyle=\sum_{i=2}^{n}\left(\sum_{k=1}^{s}u_{k}(h(x_{i},y_{1})-h(x_{i}^{\prime},y_{1}))\right)+\sum_{j=1}^{n}\left(\sum_{k=1}^{s}u_{k}(h(x_{1},y_{j})-h(x_{1}^{\prime},y_{j}^{\prime}))\right)$
$\displaystyle\leqslant
8n^{-1}\left\|u\right\|\sup_{x,y}\left\|\mathbf{C}(x,y)\right\|\left\|p\right\|_{\infty}$

So, $f(z_{1},z_{2},\ldots,z_{n})-f(z_{1}^{\prime},z_{2},\ldots,z_{n})\leqslant
8n^{-1}\left\|p\right\|_{\infty}$. and

$\left\|\Phi_{n}^{*}P-\Phi^{*}p\alpha\otimes\beta\right\|\lesssim\frac{1+m}{\sqrt{n}}$

with probability at least $1-2\exp(-2m^{2}/(64\left\|p\right\|_{\infty}^{2}))$

∎

######  Remark 22 (Adjusting for the centred cost parametrization).

In the above two propositions, we compare
$\Phi_{n}^{*}P=\sum_{i,j}\mathbf{C}_{k}(x_{i},y_{j})P_{i,j}$ with
$\Phi^{*}\pi=\int\mathbf{C}_{k}(x,y)d\pi(x,y)$ for $P$ being a discretized
version of $\pi$. The delicate issue is that in ($\mathcal{P}_{n}$), $C_{k}$
is a centralized version of
$\hat{C}_{k}\triangleq(\mathbf{C}_{k}(x_{i},y_{j}))_{i,j}$. In particular,

$C_{k}=\hat{C}_{k}-\frac{1}{n}\sum_{i=1}^{n}(\hat{C}_{k})_{i,j}-\frac{1}{n}\sum_{i=1}^{n}(\hat{C}_{k})_{i,j}+\frac{1}{n^{2}}\sum_{j=1}^{n}\sum_{i=1}^{n}(\hat{C}_{k})_{i,j}.$

Note that if $P\mathds{1}=\frac{1}{n}\mathds{1}$ and
$P^{\top}\mathds{1}=\frac{1}{n}\mathds{1}$, then

$\sum_{i,j}(C_{k})_{i,j}P_{i,j}=\sum_{i,j}(\hat{C}_{k})_{i,j}P_{i,j}-\frac{1}{n^{2}}\sum_{j=1}^{n}\sum_{k=1}^{n}(\hat{C}_{k})_{k,j}$

This last term on the RHS is neglibible because
$\Phi^{*}(\alpha\otimes\beta)=0$ by assumption of $\mathbf{C}_{k}$ being
centred: Note that

$\frac{1}{n^{2}}\sum_{j=1}^{n}\sum_{k=1}^{n}(\hat{C}_{k})_{k,j}=\hat{\Phi}_{n}^{*}(\frac{1}{n^{2}}\mathds{1}_{n\times
n}).$

Applying the above proposition,

$\left\|\frac{1}{n^{2}}\sum_{j=1}^{n}\sum_{k=1}^{n}(\hat{C}_{k})_{k,j}-\Phi^{*}(\alpha\otimes\beta)\right\|\lesssim\frac{m+\log(n)}{\sqrt{n}}$

with probability at least $1-\exp(-m)$. So, up to constants, the above two
propositions can be applied even for our centralized cost $C$.

####  B.2.2 Invertibility bound

Recall our assumption that
$M\triangleq\mathbb{E}_{(x,y)\sim\alpha\otimes\beta}[\mathbf{C}(x,y)\mathbf{C}(x,y)^{\top}]$
is invertible. In Proposition 24, we bound the deviation of
$\Phi_{n}^{*}\Phi_{n}$ from $M$ in spectral norm, and hence establish that it
is invertible with high probability.

We will make use of the following matrix Bernstein inequality.

######  Theorem 23 (Matrix Bernstein).

[Tropp et al., 2015] Let $Z_{1},\ldots,Z_{n}\in\mathbb{R}^{d\times d}$ be
independent symmetric mean-zero random matrices such that
$\left\|Z_{i}\right\|\leqslant L$ for all $i\in[n]$. Then, for all $t\geqslant
0$,

$\mathbb{E}\left\|\sum_{i}Z_{i}\right\|\leqslant\sqrt{2\sigma\log(2d)}+\frac{1}{3}L\log(2d)$

where $\sigma^{2}=\left\|\sum_{i=1}^{n}\mathbb{E}[Z_{i}^{2}]\right\|$.

######  Proposition 24.

Assume that $\log(2s)+1\leqslant n$. Then,

$\left\|\frac{1}{n^{2}}\Phi_{n}^{*}\Phi_{n}-M\right\|\lesssim\sqrt{\frac{\log(2s)}{n-1}}+\frac{m}{\sqrt{n}}$

with probability at least $1-\mathcal{O}(\exp(-m^{2}))$.

###### Proof.

Recall that $\Phi_{n}A=\sum_{k}A_{k}C_{k}$, where $C_{k}$ is the centred
version of the matrix $\hat{C}_{k}=(\mathbf{C}_{k}(x_{i},y_{j}))_{i,j}$. That
is,

$(C_{k})_{i,j}=(\hat{C}_{k})_{i,j}-\frac{1}{n}\sum_{\ell}(\hat{C}_{k})_{\ell,j}-\frac{1}{n}\sum_{m}(\hat{C}_{k})_{i,m}+\frac{1}{n^{2}}\sum_{\ell}\sum_{m}(\hat{C}_{k})_{\ell,m}.$
(30)

For simplicity, we first do the proof for
$(\Phi_{n}A)_{i,j}=\sum_{k}A_{k}\mathbf{C}_{k}(x_{i},y_{j})$, and explain at
the end how to modify the proof to account for $\Phi_{n}$ using the centred
$C_{k}$.

First,

$\displaystyle\frac{1}{n^{2}}\Phi_{n}^{*}\Phi_{n}-M$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\left(\frac{1}{n}\sum_{j=1}^{n}\mathbf{C}(x_{i},y_{j})\mathbf{C}(x_{i},y_{j})^{\top}-\int\mathbf{C}(x_{i},y)\mathbf{C}(x_{i},y)^{\top}d\beta(y)\right)$
(31)
$\displaystyle+\frac{1}{n}\sum_{i=1}^{n}\int\mathbf{C}(x_{i},y)\mathbf{C}(x_{i},y)^{\top}d\beta(y)-\int\mathbf{C}(x,y)\mathbf{C}(x,y)^{\top}d\alpha(x)d\beta(y)$
(32)

To bound the two terms in (32), let
$Z_{i}\triangleq\int\mathbf{C}(x_{i},y)\mathbf{C}(x_{i},y)^{\top}d\beta(y)-\int\mathbf{C}(x,y)\mathbf{C}(x,y)^{\top}d\alpha(x)d\beta(y)$
and observe that these are iid matrices with zero mean . By matrix bernstein
with the bounds $\left\|Z_{i}\right\|\leqslant 2$ and
$\left\|Z_{i}^{2}\right\|\leqslant 4$,

$\mathbb{E}\left\|\frac{1}{n}\sum_{i=1}^{n}Z_{i}\right\|\leqslant\frac{\sqrt{8\log(2s)}}{\sqrt{n}}+\frac{2\log(2s)}{3n}\leqslant
4\sqrt{\frac{\log(2s)}{n}}$

assuming that $\log(2s)\leqslant n$.

For the two terms in (31),

$\displaystyle\mathbb{E}\left\|\frac{1}{n}\sum_{i=1}^{n}\left(\frac{1}{n}\sum_{j=1}^{n}\mathbf{C}(x_{i},y_{j})\mathbf{C}(x_{i},y_{j})^{\top}-\int\mathbf{C}(x_{i},y)\mathbf{C}(x_{i},y)^{\top}d\beta(y)\right)\right\|$
(33)
$\displaystyle\leqslant\mathbb{E}\left\|\frac{1}{n^{2}}\sum_{i=1}^{n}\mathbf{C}(x_{i},y_{i})\mathbf{C}(x_{i},y_{i})^{\top}\right\|+\mathbb{E}\left\|\frac{1}{n}\sum_{i=1}^{n}\left(\frac{1}{n}\sum_{\begin{subarray}{c}j=1\\\
j\neq
i\end{subarray}}^{n}\mathbf{C}(x_{i},y_{j})\mathbf{C}(x_{i},y_{j})^{\top}-\int\mathbf{C}(x_{i},y)\mathbf{C}(x_{i},y)^{\top}d\beta(y)\right)\right\|$
(34)
$\displaystyle\leqslant\frac{2}{n}+\frac{n-1}{n^{2}}\sum_{i=1}^{n}\mathbb{E}\left\|\frac{1}{n-1}\sum_{\begin{subarray}{c}j=1\\\
j\neq
i\end{subarray}}^{n}\left(\mathbf{C}(x_{i},y_{j})\mathbf{C}(x_{i},y_{j})^{\top}-\int\mathbf{C}(x_{i},y)\mathbf{C}(x_{i},y)^{\top}d\beta(y)\right)\right\|$
(35)

For each $i=1,\ldots n$, let
$Y_{j}=\mathbf{C}(x_{i},y_{j})\mathbf{C}(x_{i},y_{j})-\int\mathbf{C}(x_{i},y)\mathbf{C}(x_{i},y)^{\top}d\beta(y)$
and observe that conditional on $x_{i}$, $\left\\{Y_{j}\right\\}_{j\neq i}$
are iid matrices with zero mean. The matrix Bernstein inequality applied to
$\frac{1}{n-1}\sum_{j\in[n]\setminus\left\\{i\right\\}}Y_{j}$ implies that

$\mathbb{E}\left\|\frac{1}{n-1}\sum_{j\in[n]\setminus\left\\{i\right\\}}Y_{j}\right\|\leqslant\frac{\sqrt{8\log(2s)}}{\sqrt{n-1}}+\frac{2\log(2s)}{3n-3}\leqslant
4\sqrt{\frac{\log(s)}{n}}$

assuming that $\log(2s)\leqslant n-1$.

Finally, we apply Lemma 16 to
$f(z_{1},\ldots,z_{n})=\left\|\frac{1}{n^{2}}\Phi_{n}^{*}\Phi_{n}-M\right\|$.
Let
$\Phi_{n}^{\prime}u=\sum_{k}u_{k}\mathbf{C}(x_{i}^{\prime},y_{i}^{\prime})$
with $x_{i}^{\prime}=x_{i}$ and $y_{i}^{\prime}=y_{i}$ for $i\geqslant 2$. For
each vector $u$ of unit norm,

$\displaystyle\frac{1}{n^{2}}$
$\displaystyle\langle(\Phi_{n}^{*}\Phi_{n}-(\Phi_{n}^{\prime})^{*}\Phi_{n}^{\prime})u,\,u\rangle=\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\left\lvert\mathbf{C}(x_{i},y_{j})^{\top}u\right\rvert^{2}-\left\lvert\mathbf{C}(x_{i}^{\prime},y_{j}^{\prime})^{\top}u\right\rvert^{2}$
(36)
$\displaystyle=\frac{1}{n^{2}}\sum_{i=1}^{n}\left\lvert\mathbf{C}(x_{i},y_{1})^{\top}u\right\rvert^{2}-\left\lvert\mathbf{C}(x_{i}^{\prime},y_{1}^{\prime})^{\top}u\right\rvert^{2}+\frac{1}{n^{2}}\sum_{j=2}^{n}\left\lvert\mathbf{C}(x_{1},y_{j})^{\top}u\right\rvert^{2}-\left\lvert\mathbf{C}(x_{1}^{\prime},y_{j}^{\prime})^{\top}u\right\rvert^{2}\leqslant
4n^{-1}.$ (37)

So,

$\left\|\frac{1}{n^{2}}\Phi_{n}^{*}\Phi_{n}-M\right\|\leqslant
8\sqrt{\frac{\log(2s)}{n-1}}+\frac{m}{\sqrt{n}}$

with probability at least $1-\exp(-2m^{2}/16)$.

To conclude this proof, we now discuss how to modify the above proof in the
case of the centred cost $C_{k}$, that is $\Phi_{n}A=\sum_{k}A_{k}C_{k}$ where
$C_{k}$ is as defined in (30). Note that in this case,

$\displaystyle\frac{1}{n^{2}}(\Phi_{n}^{*}\Phi_{n})_{k,\ell}$
$\displaystyle=\frac{1}{n^{2}}\sum_{i,j=1}^{n}\mathbf{C}_{k}(x_{i},y_{j})\mathbf{C}_{\ell}(x_{i},y_{j})+\frac{1}{n^{4}}\left(\sum_{p,q=1}^{n}\mathbf{C}_{\ell}(x_{p},y_{q})\right)\left(\sum_{i,j=1}^{n}\mathbf{C}_{k}(x_{i},y_{j})\right)$
(38)
$\displaystyle-\frac{1}{n^{3}}\sum_{j=1}^{n}\left(\sum_{p=1}^{n}\mathbf{C}_{\ell}(x_{p},y_{j})\right)\left(\sum_{i=1}^{n}\mathbf{C}_{k}(x_{i},y_{j})\right)-\frac{1}{n^{3}}\sum_{i=1}^{n}\left(\sum_{q=1}^{n}\mathbf{C}_{\ell}(x_{i},y_{q})\right)\left(\sum_{j=1}^{n}\mathbf{C}_{k}(x_{i},y_{j})\right)$
(39)

We already know from the previous arguments that

$\mathbb{E}\left\|\frac{1}{n^{2}}\sum_{i,j=1}^{n}\mathbf{C}_{k}(x_{i},y_{j})\mathbf{C}_{\ell}(x_{i},y_{j})-M\right\|=\mathcal{O}(\sqrt{\log(s)/n}).$

For the last term in (38), let
$\Lambda=\left\\{(i,j,p,q)\;;\;i\in[n],j\in[n]\setminus\left\\{i\right\\},\;p\in[n]\setminus\left\\{i,j\right\\},q\in[n]\setminus\left\\{i,j,p\right\\}\right\\}$.
Note that $\Lambda$ has $n(n-1)(n-2)(n-3)$ terms, and
$\Lambda^{c}=[n]\times[n]\times[n]\times[n]$ has $\mathcal{O}(n^{3})$ terms.
Therefore, we can write

$\displaystyle\mathbb{E}$
$\displaystyle\left\|\frac{1}{n^{4}}\left(\sum_{p,q=1}^{n}\mathbf{C}(x_{p},y_{q})\right)\left(\sum_{i,j=1}^{n}\mathbf{C}(x_{i},y_{j})^{\top}\right)\right\|$
$\displaystyle\leqslant\frac{1}{n^{4}}\mathbb{E}\left\|\sum_{(i,j,p,q)\in\Lambda}\mathbf{C}(x_{p},y_{q})\mathbf{C}(x_{i},y_{j})^{\top}\right\|+\mathbb{E}\left\|\frac{1}{n^{4}}\sum_{(i,j,p,q)\in\Lambda^{c}}\mathbf{C}(x_{p},y_{q})\mathbf{C}(x_{i},y_{j})^{\top}\right\|$
$\displaystyle\leqslant\frac{1}{n^{4}}\mathbb{E}\left\|\sum_{(i,j,p,q)\in\Lambda}\mathbf{C}(x_{p},y_{q})\mathbf{C}(x_{i},y_{j})^{\top}\right\|+\mathcal{O}(n^{-1})$

where the second inequality is because there are $\mathcal{O}(n^{3})$ terms in
$\Lambda^{c}$ and we used the bound that $\left\|\mathbf{C}(x,y)\right\|=1$.
Moreover,

$\displaystyle\frac{1}{n^{4}}\mathbb{E}\left\|\sum_{(i,j,p,q)\in\Lambda}\mathbf{C}(x_{p},y_{q})\mathbf{C}(x_{i},y_{j})^{\top}\right\|\leqslant\frac{n-3}{n^{4}}\sum_{i}\sum_{j\not\in\left\\{i\right\\}}\sum_{p\not\in\left\\{i,j\right\\}}\mathbb{E}\left\|\frac{1}{n-3}\sum_{q\not\in\left\\{i,j,p\right\\}}\mathbf{C}(x_{p},y_{q})\mathbf{C}(x_{i},y_{j})^{\top}\right\|.$

Note that conditional on $i,j,p$,

$\mathbb{E}[\frac{1}{n-3}\sum_{q\not\in\left\\{i,j,p\right\\}}\mathbf{C}(x_{p},y_{q})\mathbf{C}(x_{i},y_{j})^{\top}]=\int\mathbf{C}(x_{p},y)d\beta(y)\mathbf{C}(x_{i},y_{j})^{\top}=0$

by assumption that $\mathbf{C}(x,y)d\beta(y)=0$. So, we can apply Matrix
Bernstein as before to show that

$\mathbb{E}[\frac{1}{n^{4}}\mathbb{E}\left\|\sum_{(i,j,p,q)\in\Lambda}\mathbf{C}(x_{p},y_{q})\mathbf{C}(x_{i},y_{j})^{\top}\right\|]\lesssim\sqrt{\log(s)n^{-1}}.$

A similar argument can be applied to handle the two terms in (39), and so,

$\mathbb{E}\left\|\frac{1}{n^{2}}\Phi_{n}^{*}\Phi_{n}-M\right\|\lesssim\sqrt{\log(s)/n}.$

The high probability bound can now be derived using Lemma 16 as before.

∎

##  Appendix C 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
$\Sigma_{\alpha}=\mathrm{Id}$ and $\Sigma_{\beta}=\mathrm{Id}$):

$\tilde{W}(A)\triangleq\sup_{\Sigma}\langle
A,\,\Sigma\rangle+\frac{\varepsilon}{2}\log\det(\mathrm{Id}-\Sigma^{\top}\Sigma).$
(40)

To retrieve the Hessian of the original function note that, since
$\log\det(\Sigma_{\beta}-\Sigma^{\top}\Sigma_{\alpha}^{-1}\Sigma)=\log\det(\Sigma_{\beta})+\log\det(\mathrm{Id}-\Sigma_{\beta}^{-\frac{1}{2}}\Sigma^{\top}\Sigma_{\alpha}^{-1}\Sigma\Sigma_{\beta}^{-\frac{1}{2}})$,
a change of variable
$\tilde{\Sigma}\triangleq\Sigma_{\alpha}^{-\frac{1}{2}}\Sigma\Sigma_{\beta}^{-\frac{1}{2}}$,
shows that
$W(A)=\tilde{W}(\Sigma_{\alpha}^{\frac{1}{2}}A\Sigma_{\beta}^{\frac{1}{2}})$
and, hence,

$\nabla^{2}W(A)=(\Sigma_{\beta}^{\frac{1}{2}}\otimes\Sigma_{\alpha}^{\frac{1}{2}})\nabla^{2}\tilde{W}(\Sigma_{\alpha}^{\frac{1}{2}}A\Sigma_{\beta}^{\frac{1}{2}})(\Sigma_{\beta}^{\frac{1}{2}}\otimes\Sigma_{\alpha}^{\frac{1}{2}}).$
(41)

By the envelope theorem, $\nabla\tilde{W}(A)=\Sigma$, where $\Sigma$ is the
maximizer of (40), thus reducing the computation of $\nabla^{2}\tilde{W}(A)$
to differentiating the optimality condition of $\Sigma$, _i.e.,_

$A=\varepsilon^{-1}\Sigma(\mathrm{Id}-\Sigma^{\top}\Sigma)^{-1}.$ (42)

Recall also that such a $\Sigma$ has an explicit formula given in (4).

######  Lemma 25.

$\displaystyle\nabla^{2}\tilde{W}(A)$
$\displaystyle=\varepsilon\left(\varepsilon^{2}(\mathrm{Id}-\Sigma^{\top}\Sigma)^{-1}\otimes(\mathrm{Id}-\Sigma\Sigma^{\top})^{-1}+(A^{\top}\otimes
A)\mathbb{T}\right)^{-1}$ (43)

and $\Sigma$ is as in (42) and $\mathbb{T}$ is the linear map defined by
$\mathbb{T}\mbox{vec}(X)=\mbox{vec}(X^{\top})$.

###### Proof of Lemma 25.

Define
$G(\Sigma,A)\triangleq\Sigma(\mathrm{Id}-\Sigma^{\top}\Sigma)^{-1}-\varepsilon^{-1}A$.
Note that a maximizer $\Sigma$ of (40) satisfies $G(\Sigma,A)=0$. Moreover,
$\partial_{\Sigma}G$ is invertible at such $(\Sigma,A)$ because this is the
Hessian of $\log\det(\mathrm{Id}-\Sigma^{\top}\Sigma)$, a strictly concave
function. By the implicit function theorem, there exists a function
$f:A\mapsto\Sigma$ such that $G(f(A),A)=0$ and

$\nabla
f=(\partial_{\Sigma}G)^{-1}\partial_{A}G=\varepsilon^{-1}(\partial_{\Sigma}G)^{-1}.$

It remains to compute $\partial_{\Sigma}G$ at $(f(A),A)$:

At an arbitrary point $(\Sigma,A)$ we have

$\displaystyle\partial_{\Sigma}G$
$\displaystyle=(\mathrm{Id}-\Sigma^{\top}\Sigma)^{-1}\otimes\mathrm{Id}+(\mathrm{Id}\otimes\Sigma)\partial_{\Sigma}(\mathrm{Id}-\Sigma^{\top}\Sigma)^{-1}$
$\displaystyle=(\mathrm{Id}-\Sigma^{\top}\Sigma)^{-1}\otimes\mathrm{Id}+(\mathrm{Id}\otimes\Sigma)\left((\mathrm{Id}-\Sigma^{\top}\Sigma)^{-1}\otimes(\mathrm{Id}-\Sigma^{\top}\Sigma)^{-1}\right)\left((\Sigma^{\top}\otimes\mathrm{Id})\mathbb{T}+(\mathrm{Id}\otimes\Sigma^{\top})\right)$

and at $(f(A),A)$, _i.e,_ with
$\Sigma(\mathrm{Id}-\Sigma^{\top}\Sigma)^{-1}=\varepsilon^{-1}A$, we can
further simplify to

$\partial_{\Sigma}G=(\mathrm{Id}-\Sigma^{\top}\Sigma)^{-1}\otimes(\mathrm{Id}+\Sigma(\mathrm{Id}-\Sigma^{\top}\Sigma)^{-1}\Sigma^{\top}))+\varepsilon^{-2}(A^{\top}\otimes
A)\mathbb{T}$

By the woodbury matrix formula,

$(\mathrm{Id}+\Sigma(\mathrm{Id}-\Sigma^{\top}\Sigma)^{-1}\Sigma^{\top})=(\mathrm{Id}-\Sigma\Sigma^{\top})^{-1}=\mathrm{Id}+A\Sigma^{\top}.$

So,

$\partial_{\Sigma}G=(\mathrm{Id}-\Sigma^{\top}\Sigma)^{-1}\otimes(\mathrm{Id}-\Sigma\Sigma^{\top})^{-1}+\varepsilon^{-2}(A^{\top}\otimes
A)\mathbb{T},$

thus concluding the proof. ∎

We remark that, from the connection between $W(A)$ and $\tilde{W}(A)$, _i.e.,_
41, we obtain Lemma 8 as a corollary.

###  C.1 Limit cases

##### SVD representation of the covariance

To derive the limiting expressions, we make an observation on the singular
value decomposition of $\Sigma$: Let the singular value decomposition of $A$
be $A=UDV^{\top}$, where $D$ is the diagonal matrix with positive entries
$d_{i}$. Note that
$\Delta=(\mathrm{Id}+\frac{\varepsilon^{2}}{4}(A^{\top}A)^{\dagger})^{\frac{1}{2}}=V(\mathrm{Id}+\frac{\varepsilon^{2}}{4}D^{-2})^{\frac{1}{2}}V^{\top}$.
Moreover, $\Delta$ and $A^{\top}A$ commute, so

$\displaystyle\Sigma$
$\displaystyle=A\Delta((\Delta^{2}A^{\top}A)^{\dagger})^{\frac{1}{2}}\Delta-\frac{\varepsilon}{2}A^{T,\dagger}=U\left(\left(\mathrm{Id}+\frac{\varepsilon^{2}}{4}D^{\dagger,2}\right)^{\frac{1}{2}}-\frac{\varepsilon}{2}D^{\dagger}\right)V^{\top}=U\tilde{D}V^{\top},$
(44)

where $\tilde{D}$ is the diagonal matrix with diagonal entries

$\tilde{d}_{i}=\sqrt{\left(1+\frac{\varepsilon^{2}}{4d_{i}^{2}}\right)}-\frac{\varepsilon}{2}\frac{1}{d_{i}}.$
(45)

####  C.1.1 Link with Lasso: Limit as $\varepsilon\to\infty$

Note that

$\displaystyle\tilde{d}_{i}$
$\displaystyle=\frac{\varepsilon}{2d_{i}}\sqrt{1+\frac{4d_{i}^{2}}{\varepsilon^{2}}}-\frac{\varepsilon}{2d_{i}}=\frac{d_{i}}{\varepsilon}-\frac{d_{i}^{3}}{\varepsilon^{3}}+\mathcal{O}(\varepsilon^{-5})\to
0,\qquad\varepsilon\to\infty.$

It follows that $\lim_{\varepsilon\to\infty}\Sigma=0$ and hence,
$\varepsilon\nabla^{2}\tilde{W}(A)\to\mathrm{Id}$. So, the certificate
converges to

$(\Sigma_{\beta}\otimes\Sigma_{\alpha})_{(:,I)}\big{(}(\Sigma_{\beta}\otimes\Sigma_{\alpha})_{(I,I)}\big{)}^{-1}\operatorname{sign}({\widehat{A}})_{I}.$
(46)

###### Proof of Proposition 9.

The iOT problem approaches a Lasso problem as $\varepsilon\to\infty$. Recall
that in the Gaussian setting, the iOT problem is of the form

$\operatorname*{arg\,min}_{A}\mathcal{F}(A)=\lambda\left\|A\right\|_{1}+\langle
A,\,\Sigma_{\varepsilon,A}-\Sigma_{\varepsilon,{\widehat{A}}}\rangle+\frac{\varepsilon}{2}\log\det\left(\Sigma_{\beta}-\Sigma_{\varepsilon,A}^{\top}\Sigma_{\alpha}^{-1}\Sigma_{\varepsilon,A}\right)$
(47)

where $\Sigma_{\varepsilon,A}$ satisfies
$\Sigma_{\beta}-\Sigma_{\varepsilon,A}^{\top}\Sigma_{\alpha}^{-1}\Sigma_{\varepsilon,A}=\varepsilon\Sigma_{\alpha}^{-1}\Sigma_{\varepsilon,A}A^{-1}.$
So, we can write

$\operatorname*{arg\,min}_{A}\mathcal{F}(A)\triangleq\lambda\left\|A\right\|_{1}+\langle
A,\,\Sigma_{\varepsilon,A}-\Sigma_{\varepsilon,{\widehat{A}}}\rangle+\frac{\varepsilon}{2}\log\det\left(\varepsilon\Sigma_{\alpha}^{-1}\Sigma_{\varepsilon,A}A^{-1}\right)$
(48)

Let

$X\triangleq\Sigma_{\alpha}^{\frac{1}{2}}A\Sigma_{\beta}^{\frac{1}{2}}\quad\text{and}\quad\tilde{\Sigma}_{\varepsilon,A}=\Sigma_{\alpha}^{-\frac{1}{2}}\Sigma_{\varepsilon,A}\Sigma_{\beta}^{-\frac{1}{2}}.$

From (45), if $X$ has svd decomposition
$W=U\operatorname{diag}(d_{i})V^{\top}$, then

$\tilde{\Sigma}_{\varepsilon,A}=U\tilde{D}V^{\top},\quad\text{where}\quad\tilde{D}=\operatorname{diag}(\tilde{d}_{i})$

and

$\displaystyle\tilde{d}_{i}=\frac{d_{i}}{\varepsilon}-\frac{d_{i}^{3}}{\varepsilon^{3}}+\mathcal{O}(\varepsilon^{-5}).$

So,

$\displaystyle\log\det(\varepsilon\Sigma_{\alpha}^{-1}\Sigma_{\varepsilon,A}A^{-1})$
$\displaystyle=\log\det\left(\varepsilon\Sigma_{\alpha}^{-\frac{1}{2}}\tilde{\Sigma}_{\varepsilon,A}X^{-1}\Sigma_{\alpha}^{\frac{1}{2}}\right)$
$\displaystyle=\log\det\left(U\operatorname{diag}\left(1-d_{i}^{2}/\varepsilon^{2}+\mathcal{O}(\varepsilon^{-4})\right)U^{\top}\right)$
$\displaystyle=\log\det\left(\operatorname{diag}\left(1-d_{i}^{2}/\varepsilon^{2}+\mathcal{O}(\varepsilon^{-4})\right)\right)=-\frac{1}{\varepsilon^{2}}\left\|X\right\|_{F}^{2}+\mathcal{O}(\varepsilon^{-4}).$

Also,

$\varepsilon\langle\Sigma_{\varepsilon,A}-\Sigma_{\varepsilon,{\widehat{A}}},\,A\rangle=\varepsilon\langle\tilde{\Sigma}_{\varepsilon,A}-\tilde{\Sigma}_{\varepsilon,{\widehat{A}}},\,X\rangle=\langle
X-X_{0},\,X\rangle+\mathcal{O}(\varepsilon^{-2}).$

So, assuming that $\lambda=\lambda_{0}/\varepsilon$,

$\displaystyle\varepsilon\mathcal{F}(A)$
$\displaystyle=\lambda_{0}\left\|A\right\|_{1}+\langle
X-X_{0},\,X\rangle-\left\|X\right\|_{F}^{2}+\mathcal{O}(\varepsilon^{-2})$
(49)
$\displaystyle=\lambda_{0}\left\|A\right\|_{1}+\left\|(\Sigma_{\beta}^{\frac{1}{2}}\otimes\Sigma_{\alpha}^{\frac{1}{2}})(A-{\widehat{A}})\right\|_{F}^{2}-\frac{1}{2}\left\|(\Sigma_{\beta}^{\frac{1}{2}}\otimes\Sigma_{\alpha}^{\frac{1}{2}}){\widehat{A}}\right\|^{2}+\mathcal{O}(\varepsilon^{-2})$
(50)

The final statement on convergence of minimizers follows by Gamma convergence.
∎

####  C.1.2 Link with Graphical Lasso

In the special case where the covariances are the identity
($\Sigma_{\alpha}=\mathrm{Id}$ and $\Sigma_{\beta}=\mathrm{Id}$) and $A$ is
symmetric positive definite we have that $\Sigma$ is also positive definite
and Galichon’s formula (4) holds (since $A$ is invertible) and hence the
Hessian reduces to

$\displaystyle\left(\frac{1}{\varepsilon}\nabla^{2}\tilde{W}(A)\right)^{-1}=(A\otimes
A)\left(\Sigma^{-1}\otimes\Sigma^{-1}+\mathbb{T}\right).$ (51)

Moreover, if $A$ admits an eigenvalue decomposition $A=UDU^{\top}$, then
$\Sigma$ admits an eigenvalue decomposition $\Sigma=U\tilde{D}U^{\top}$ with
entries of $\tilde{D}$ given by (45). Note that it follows that
$\lim_{\varepsilon\to 0}\Sigma=\mathrm{Id}$ and, hence, $\lim_{\varepsilon\to
0}\left(\Sigma^{-1}\otimes\Sigma^{-1}+\mathbb{T}\right)=\mathrm{Id}+\mathbb{T}$,
a singular matrix with the kernel being the set of asymmetric matrices. So,
the limit does not necessarily exist as $\varepsilon\to 0$. However, in the
special case where $A$ is _symmetric positive definite_ , one can show that
the certificates remain well defined as $\varepsilon\to 0$: Let $S_{+}$ be the
set of matrices $\psi$ such that $\psi$ is symmetric and $S_{-}$ be the set of
matrices $\psi$ such that $\psi$ is anti-symmetric. Then,
$(\nabla^{2}\tilde{W}(A))^{-1}(S_{+})\subset S_{+}$ and
$(\nabla^{2}\tilde{W}(A))^{-1}(S_{-})\subset S_{-}$. Moreover,

$\displaystyle\left(\frac{1}{\varepsilon}\nabla^{2}\tilde{W}(A)\right)^{-1}\restriction_{S_{+}}$
$\displaystyle=(A\otimes
A)\left(\Sigma^{-1}\otimes\Sigma^{-1}+\mathrm{Id}\right),$
$\displaystyle\left(\frac{1}{\varepsilon}\nabla^{2}\tilde{W}(A)\right)^{-1}\restriction_{S_{+}}$
$\displaystyle=(A\otimes
A)\left(\Sigma^{-1}\otimes\Sigma^{-1}-\mathrm{Id}\right).$

Since the symmetry of $A$ implies the symmetry of $\mbox{sign}(A)$, we can
replace the Hessian given in (51) by $(A\otimes
A)(\Sigma^{-1}\otimes\Sigma^{-1}+\mathrm{Id})$ and, hence, since
$\lim_{\varepsilon\to 0}\Sigma=\mathrm{Id}$, the limit as $\varepsilon\to 0$
is

$\lim_{\varepsilon\to 0}z_{\varepsilon}=(A^{-1}\otimes
A^{-1})_{(:,I)}\big{(}(A^{-1}\otimes
A^{-1})_{(I,I)}\big{)}^{-1}\operatorname{sign}(A)_{I}.$ (52)

This coincides precisely with the certificate of the graphical lasso:

$\operatorname*{arg\,min}_{\Theta\succeq 0}\langle
S,\,\Theta\rangle-\log\det(\Theta)+\lambda\left\|\Theta\right\|_{1}.$

###### Proof of Proposition 10.

The iOT problem with identity covariances restricted to the set of positive
semi-definite matrices has the form

$\operatorname*{arg\,min}_{A\succeq
0}\mathcal{F}(A)_{\varepsilon,\lambda}(A),\quad\text{where}\quad\mathcal{F}(A)_{\varepsilon,\lambda}(A)\triangleq\lambda\left\|A\right\|_{1}+\langle
A,\,\Sigma_{\varepsilon,A}-\hat{\Sigma}\rangle+\frac{\varepsilon}{2}\log\det(\mathrm{Id}-\Sigma_{\varepsilon,A}^{\top}\Sigma_{\varepsilon,A}),$
(53)

where
$I-\Sigma_{\varepsilon,A}^{\top}\Sigma_{\varepsilon,A}=\varepsilon\Sigma_{\varepsilon,A}A^{-1}.$
From the singular value decomposition of $\Sigma_{\varepsilon,A}$,_i.e.,_
(44), we see that if $A$ is symmetric positive definite, then so is
$\Sigma_{\varepsilon,A}$. Plugging the optimality condition, _i.e.,_
$I-\Sigma_{\varepsilon,A}^{\top}\Sigma_{\varepsilon,A}=\varepsilon\Sigma_{\varepsilon,A}A^{-1}$,
into (53), we obtain

$\displaystyle\operatorname*{arg\,min}_{A\succeq
0}\lambda\left\|A\right\|_{1}+\langle
A,\,\Sigma_{\varepsilon,A}-\hat{\Sigma}\rangle+\frac{\varepsilon}{2}\log\det(\varepsilon
A^{-1}\Sigma_{\varepsilon,A})$ $\displaystyle=$
$\displaystyle\operatorname*{arg\,min}_{A\succeq
0}\lambda\left\|A\right\|_{1}-\frac{\varepsilon}{2}\log\det(A/\varepsilon)+\varepsilon\langle
A,\,(\Sigma_{\varepsilon,A}-\hat{\Sigma})/\varepsilon\rangle+\frac{\varepsilon}{2}\log\det(\Sigma_{\varepsilon,A})$

Let $\lambda=\varepsilon\lambda_{0}$ for some $\lambda_{0}>0$. Then, removing
the constant $\frac{\varepsilon}{2}\log\det(\varepsilon\mathrm{Id})$ term and
factoring out $\varepsilon$, the problem is equivalent to

$\operatorname*{arg\,min}_{A\succeq
0}\lambda_{0}\left\|A\right\|_{1}-\frac{1}{2}\log\det(A)+\langle
A,\,\varepsilon^{-1}(\Sigma_{\varepsilon,A}-\hat{\Sigma})\rangle+\frac{1}{2}\log\det(\Sigma_{\varepsilon,A})$

Assume that $\hat{\Sigma}=\Sigma_{\varepsilon,{\widehat{A}}}$. From the
expression for the singular values of $\Sigma_{\varepsilon,A}$ in (45), note
that

$\Sigma_{\varepsilon,A}=\mathrm{Id}-\frac{\varepsilon}{2}A^{-1}+\mathcal{O}(\varepsilon^{2}).$

So, $\lim_{\varepsilon\to
0}(\Sigma_{\varepsilon,A}-\Sigma_{\varepsilon,{\widehat{A}}})/\varepsilon=-\frac{1}{2}\left(A^{-1}-{\widehat{A}}^{-1}\right)$.
The objective converges pointwise to

$\lambda_{0}\left\|A\right\|_{1}-\frac{1}{2}\log\det(A)+\frac{1}{2}\langle
A,\,{\widehat{A}}^{-1}-A^{-1}\rangle,$

and the statement is then a direct consequence of Gamma convergence. ∎

##  Appendix D Large scale $\ell^{1}$-iOT solver

Recall that the iOT optimization problem, recast over the dual potentials for
empirical measures, reads

$\inf_{A,F,G}\frac{1}{n}\sum_{i}F_{i}+G_{i}+(\Phi_{n}A)_{i,i}+\frac{\varepsilon}{2n^{2}}\sum_{i,j}\exp\left(\frac{2}{\varepsilon}(F_{i}+G_{j}+(\Phi_{n}A)_{i,j}\right)+\lambda\left\|A\right\|_{1}.$

In order to obtain a better-conditioned optimization problem, in line with
[Cuturi and 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}\frac{1}{n}\sum_{i}F_{i}+(\Phi_{n}A)_{i,i}+\frac{\varepsilon}{n}\sum_{i}\log\frac{1}{n}\sum_{j}\exp\left(\frac{2}{\varepsilon}(F_{i}+(\Phi_{n}A)_{i,j}\right)+\lambda\left\|A\right\|_{1}.$

Following [Poon and Peyré, 2021], which proposes a state of the art Lasso
solver, the last step is to use the following Hadamard product over-
parameterization of the $\ell^{1}$ norm

$\left\|A\right\|_{1}=\min_{U\odot
V}\frac{\left\|U\right\|^{2}_{2}}{2}+\frac{\left\|V\right\|^{2}_{2}}{2}.$

where the Hadamard product is $U\odot V\triangleq(U_{i}V_{i})_{i}$, to obtain
the final optimization problem

$\inf_{A,U,V}\frac{1}{n}\sum_{i}F_{i}+(\Phi_{n}(U\odot
V))_{i,i}+\frac{\varepsilon}{n}\sum_{i}\log\frac{1}{n}\sum_{j}\exp\left(\frac{2}{\varepsilon}(F_{i}+(\Phi_{n}(U\odot
V))_{i,j}\right)+\frac{\lambda}{2}\left\|U\right\|^{2}_{2}+\frac{\lambda}{2}\left\|V\right\|^{2}_{2}.$

This is a smooth optimization problem, for which we employ a quasi-Newton
solver (L-BFGS). Although it is non-convex, as demonstrated in [Poon and
Peyré, 2021], the non-convexity is benign, ensuring the solver always
converges to a global minimizer, $(F^{\star},U^{\star},V^{\star})$, of the
functional. From this, one can reconstruct the cost parameter,
$A^{\star}\triangleq U^{\star}\odot V^{\star}$.