# Outlier-Robust Optimal Transport

Debarghya Mukherjee  Aritra Guha  Justin Solomon  Yuekai Sun  Mikhail
Yurochkin

###### Abstract

Optimal transport (OT) measures distances between distributions in a way that
depends on the geometry of the sample space. In light of recent advances in
computational OT, OT distances are widely used as loss functions in machine
learning. Despite their prevalence and advantages, OT loss functions can be
extremely sensitive to outliers. In fact, a single adversarially-picked
outlier can increase the standard $W_{2}$-distance arbitrarily. To address
this issue, we propose an outlier-robust formulation of OT. Our formulation is
convex but challenging to scale at a first glance. Our main contribution is
deriving an _equivalent_ formulation based on cost truncation that is easy to
incorporate into modern algorithms for computational OT. We demonstrate the
benefits of our formulation in mean estimation problems under the Huber
contamination model in simulations and outlier detection tasks on real data.

Machine Learning, ICML

Appendix

Appendix

##  1 Introduction

Optimal transport (OT) is a fundamental problem in applied mathematics. In its
original form (Monge, 1781), the problem seeks the minimum-cost way to
transport mass from a probability distribution $\mu$ on $\mathcal{X}$ to
another distribution $\nu$ on $\mathcal{X}$. In its original form, Monge’s
problem proved hard to study, and Kantorovich (1942) relaxed Monge’s
formulation of the optimal transport problem to

$\text{OT}(\mu,\nu)\triangleq\min_{\Pi\in\mathcal{F}(\mu,\nu)}{\mathbb{E}}_{(X_{1},X_{2})\sim\Pi}\big{[}c(X_{1},X_{2})\big{]},$
(1.1)

where $\mathcal{F}(\mu,\nu)$ is the set of couplings between $\mu$ and $\nu$
(probability distributions on $\mathcal{X}\times\mathcal{X}$ whose marginals
are $\mu$ and $\nu$) and $c$ is a cost function. In this paper, we assume
$c(x,y)\geq 0$ and $c(x,x)=0$. Compared to other common measure of distance
between probability distributions (e.g.$d$-divergences), optimal transport
uniquely depends on the geometry of the sample space (through the cost
function).

Recent advancements in optimization for optimal transport (Cuturi, 2013;
Solomon et al., 2015; Genevay et al., 2016; Seguy et al., 2018) enable its
broad adaptation in machine learning applications where geometry of the data
is important; see (Peyré & Cuturi, 2018) for a survey. Optimal transport has
found applications in natural language processing (Kusner et al., 2015; Huang
et al., 2016; Alvarez-Melis & Jaakkola, 2018; Yurochkin et al., 2019),
generative modeling (Arjovsky et al., 2017), clustering (Ho et al., 2017),
domain adaptation (Courty et al., 2014, 2017), large-scale Bayesian modeling
(Srivastava et al., 2018), anomaly detection (Tong et al., 2020), and many
other domains.

Many applications use OT as a loss in an optimization problem of the form:

$\textstyle\theta\in\operatorname*{arg\,min}_{\theta\in\Theta}\text{OT}(\mu_{n},\nu_{\theta}),$
(1.2)

where $\\{\nu_{\theta}\\}_{\theta\in\Theta}$ is a collection of parametric
models and $\mu_{n}$ is the empirical distribution of the samples. Such
estimators are called _minimum Kantorovich estimators (MKE)_ (Bassetti et al.,
2006). They are popular alternatives to likelihood-based estimators,
especially in generative modeling. For example, when $\text{OT}(\cdot,\cdot)$
is the Wasserstein-1 distance and $\nu_{\theta}$ is a generator parameterized
by a neural network with weights $\theta$, equation 1.2 corresponds to the
Wasserstein GAN (Arjovsky et al., 2017).

One drawback of optimal transport is its sensitivity to outliers. Because
_all_ the mass in $\mu$ must be transported to $\nu$, a small fraction of
outliers can have an outsized impact. For statistics and machine learning
applications in which the data is corrupted or noisy, this is a major issue.
For example, the poor performance of Wasserstein GANs in the presence of
outliers was noted in the recent works on outlier-robust generative learning
with $f$-divergence GANs (Chao et al., 2018; Wu et al., 2020). The problem of
outlier-robustness in MKE has not been studied except in two recent works
proposing changes to the OT formulation that are challenging to handle
computationally (Staerman et al., 2020; Balaji et al., 2020). Our goal is to
derive an outlier-robust OT formulation compatible with existing efficient
computational OT methods (Peyré & Cuturi, 2018).

In this paper, we propose a modification of OT to address its sensitivity to
outliers. Our formulation can be used as a loss in equation 1.2, so that it is
robust to a small fraction of outliers in the data. For simplicity, we
consider the ${\epsilon}$-contamination model (Huber & Ronchetti, 2009). Let
$\nu_{\theta_{0}}$ be a member of a parametric family
$\\{\nu_{\theta}:\theta\in\Theta\\}$ and let

$\mu=(1-{\epsilon})\nu_{\theta_{0}}+{\epsilon}\tilde{\nu},$

where $\mu$ is the data-generating distribution, ${\epsilon}>0$ is the
fraction of outliers, and $\tilde{\nu}$ is the distribution of the outliers.
Although the fraction of outliers is capped at ${\epsilon}$, the value of the
outliers can be arbitrary, so their effect on the optimal transport problem
can be arbitrarily large. We modify the problem so that it is more robust to
such outliers, targeting the downstream application of learning $\theta_{0}$
from (samples from) $\mu$ in the ${\epsilon}$-contamination model.

Our main contributions are as follows:

  * •

We propose a robust OT formulation suitable for statistical estimation in the
${\epsilon}$-contamination model using MKE.

  * •

We show that our formulation is equivalent to the original OT problem with a
clipped transport cost. This connection enables us to leverage the voluminous
literature on computational optimal transport to develop efficient algorithm
to perform MKE robust to outliers.

  * •

Our formulation enables a new application of optimal transport: outlier
detection in data.

##  2 Problem Formulation

###  2.1 Robust OT for MKE

To promote outlier-robustness in MKE, we need to allow the corresponding OT
problem to ignore outliers in the data distribution $\mu$. The
${\epsilon}$-contamination model imposes a cap on the fraction of outliers, so
it is not hard to see that $\|\mu-\nu_{\theta_{0}}\|_{\
\mathrm{TV}}\leq{\epsilon}$, where $\|\cdot\|_{\ \mathrm{TV}}$ is the total-
variation norm defined as
$\|\mu\|_{\text{TV}}=\int\frac{1}{2}|\mu(\mathrm{d}x)|$. This suggests we
solve a TV-constrained/regularized version of equation 1.2:

$\displaystyle\min_{\theta\in\Theta,\tilde{\mu}}$
$\displaystyle\text{OT}(\tilde{\mu},\nu_{\theta})$ (2.1) $\displaystyle{\sf
subject\ to}$ $\displaystyle\|\mu-\tilde{\mu}\|_{\
\mathrm{TV}}\leq{\epsilon};\ \ \tilde{\mu}\in\mathbf{P},$

where $\mathbf{P}$ is the set of all probability measures. The constrained
version, however, suffers from identification issues. It cannot distinguish
between “clean” distributions within TV distance ${\epsilon}$ of
$\nu_{\theta_{0}}$. To see this, fix $\theta\in\Theta$, and note that the
optimal value of equation 2.1 is zero whenever $\mu$ is within ${\epsilon}$-TV
distance of $\nu_{\theta}$. Thus equation 2.1 cannot distinguish between two
parameter values $\theta_{1}$ and $\theta_{2}$ such that
$\|\nu_{\theta_{1}}-\nu_{\theta_{2}}\|_{{\sf TV}}\leq{\epsilon}$. This makes
it unsuitable as a loss function for statistical estimation, because it cannot
lead to a consistent estimator. As an alternative, its regularized counterpart
does not suffer from this issue:

$\min_{\begin{subarray}{c}\theta\in\Theta\\\
s:\mu+s\in\mathbf{P}\end{subarray}}\text{OT}(\mu+s,\nu_{\theta})+\lambda\|s\|_{\
\mathrm{TV}},$ (2.2)

where $\lambda>0$ is a regularization parameter. Note that, the constrained
and Lagrangian formulations are equivalent, but the equivalence depends on the
two distributions in the arguments of the Wasserstein distance. In other
words, as we vary the distributions (keeping ${\epsilon}$ fixed), the
equivalent $\lambda$ will change. In our formulation, we are fixing $\lambda$
and varying the distributions, so the solution paths are different, making the
parameter identifiable. In the rest of this paper, we work with the TV-
regularized formulation equation 2.2.

The main idea of our formulation is to allow for modifications of $\mu$, while
penalizing their magnitude and ensuring that the modified $\mu$ is still a
probability measure. Below we formulate this intuition in an optimization
problem titled ROBOT (ROBust Optimal Transport):

Formulation 1:

$\displaystyle\begin{aligned} &\min_{\pi,s}&&\iint c(x,y)\ \pi(x,y)\ dx\
dy+\lambda\|s\|_{\mathrm{TV}}\\\ &{\sf subject\ to}&&\int\pi(x,y)\
dy=\mu(x)+s(x)\geq 0\\\ &&&\int\pi(x,y)\ dy=\nu(y)\\\ &&&\int
s(\mathrm{d}x)=0.\end{aligned}$ (2.3)

where $\pi$ is a density function on $\mathcal{X}\times\mathcal{X}$. The first
and the last constraints ensure that $\mu+s$ is a valid probability measure,
while $\lambda\|s\|_{\text{TV}}$ penalizes the amount of modifications in
$\mu$. We can identify exact locations of outliers in $\mu$ by inspecting
$\mu+s$, i.e. if $\mu(x)+s(x)=0$, then $x$ got eliminated and is an outlier.
We will use this property to propose an outlier detection method.

ROBOT, unlike classical OT, guarantees that an adversarially-picked outliers
cannot increase the (robust) transport distance arbitrarily. Let
$\tilde{\mu}=(1-{\epsilon})\mu+{\epsilon}\mu_{c}$, i.e., $\tilde{\mu}$ is
$\mu$ contaminated with outliers from $\mu_{c}$, and let $\nu$ be an arbitrary
measure; in MKE, $\tilde{\mu}$ is the contaminated data and $\nu$ is the model
we learn. The adversary can arbitrarily increase $\text{OT}(\tilde{\mu},\nu)$
by manipulating the outlier distribution $\mu_{c}$. For ROBOT, we have the
following bound:

######  Theorem 2.1.

Let $\tilde{\mu}=(1-{\epsilon})\mu+{\epsilon}\mu_{c}\,$ for some
${\epsilon}\in[0,1)$. Then,

$\displaystyle\mathrm{ROBOT}(\tilde{\mu},\nu)$
$\displaystyle\leq\min\left\\{\mathrm{OT}(\mu,\nu)+\lambda{\epsilon}\|\mu-\mu_{c}\|_{\mathrm{TV}},\right.$
$\displaystyle\qquad\qquad\left.\lambda\|\tilde{\mu}-\nu\|_{\mathrm{TV}},\mathrm{OT}(\tilde{\mu},\nu)\right\\},$
(2.4)

This bound has two key takeaways: since TV norm of distributions is bounded by
1, the adversary can not increase ROBOT$(\tilde{\mu},\nu)$ arbitrarily; in the
absence of outliers, ROBOT is bounded by classical OT. See Appendix B for the
proof.

#### Related work.

We note a connection between equation 2.3 and unbalanced OT (UOT) (Chizat.,
2017; Chizat et al., 2018). UOT is typically formulated by replacing the TV
norm with KL$(\mu+s|\mu)$ and adding an analogous term for $\nu$. Chizat et
al. (2018) studied entropy regularized UOT with various divergences penalizing
marginal violations. Optimization problems similar to equation 2.3 have also
been considered outside of ML (Piccoli & Rossi, 2014; Liero et al., 2018).
Balaji et al. (2020) use UOT with $\chi^{2}$-divergence penalty on marginal
violations to achieve outlier-robustness in generative modeling. Another
relevant variation of OT is partial OT (Figalli, 2010; Caffarelli & McCann,
2010). It may also be considered for outlier-robustness but has a drawback of
forcing mass destruction rather than adjusting marginals to ignore outliers
when they are present. Staerman et al. (2020) take a different path: they
replace the expectation in the Wasserstein-1 dual with a median-of-means to
promote robustness. It is unclear what is the corresponding primal problem,
making their formulation hard to interpret as an optimal transport problem.

A major challenge with the aforementioned methods, including our Formulation
1, is the large scale implementation of the optimization problem. Chizat et
al. (2018) propose a Sinkhorn-like algorithm for entropy regularized UOT, but
it is not amenable to stochastic optimization. Balaji et al. (2020) propose a
stochastic optimization algorithm based on the UOT dual, but it requires two
additional neural networks (total of four including dual potentials) to
parameterize modified marginal distributions (i.e., $\mu+s$ and analogous one
for $\nu$). Optimizing with a median-of-means in the objective function
(Staerman et al., 2020) is also challenging. The key contribution of our work
is a formulation _equivalent_ to equation 2.3, which is _easily compatible_
with the large body of classical OT optimization techniques (Cuturi, 2013;
Solomon et al., 2015; Genevay et al., 2016; Seguy et al., 2018).

#### More efficient equivalent formulation.

At a first glance, there are two issues with equation 2.3: it appears
asymmetric, and it is unclear if it can be optimized efficiently. Below we
present an _equivalent_ formulation that is free of these issues:
Formulation 2:

$\displaystyle\min_{\Pi\in\mathcal{F}^{+}(\mathbb{R}^{d}\times\mathbb{R}^{d})}$
$\displaystyle\mathbb{E}_{(X,Y)\sim\Pi}\left[C_{\lambda}(X,Y)\right]$ (2.5)
$\displaystyle{\sf subject\ to}$ $\displaystyle X\sim\mu,\ \ Y\sim\nu\,.$

where $C_{\lambda}$ is the _truncated cost_ function defined as
$C_{\lambda}(x,y)=\min\left\\{c(x,y),2\lambda\right\\}$. Looking at equation
2.5, it is not apparent that it adds robustness to MKE, but it is symmetric,
easy to combine with entropic regularization by simply truncating the cost,
and benefits from stochastic optimization algorithms (Genevay et al., 2016;
Seguy et al., 2018).

This formulation also has a distant relation to the idea of loss truncation
for achieving robustness (Shen & Sanghavi, 2019). Pele & Werman (2009)
consider the Earth Mover Distance (discrete OT) with truncated cost to achieve
computational improvements; they also mentioned its potential to promote
robustness against outlier noise but did not explore this direction.

In Section 3, we establish an _equivalence_ between the two ROBOT
formulations, equation 2.3 and equation 2.5. This equivalence allows us to
obtain an efficient algorithm based on equation 2.5 for robust MKE. We also
provide a simple procedure for computing the optimal $s$ in equation 2.3 from
the solution of equation 2.5, enabling a new OT application: outlier
detection. We verify the effectiveness of robust MKE and outlier detection in
our experiments in Section 4. Before presenting the equivalence proof,
however, we formulate the discrete analogs of the two ROBOT formulations for
their practical value.

###  2.2 Discrete ROBOT formulations

In practice, we typically encounter samples from distributions, rather then
the distributions themselves. Sampling is also built into stochastic
optimization. In this subsection, we present the discrete versions of the
ROBOT formulations. The key detail is that, in equation 2.3, $\mu,\nu$ and $s$
are all supported on $\mathbb{R}^{d}$, while in the discrete case the
empirical measures $\mu_{n}\in\Delta^{n-1}$ and $\nu_{m}\in\Delta^{m-1}$ are
supported on a set of points ($\Delta^{r}$ is the unit probability simplex in
${\mathbb{R}}^{r}$). As a result, to formulate a discrete version of equation
2.3, we need to augment $\mu_{n}$ and $\nu_{m}$ with each others’ supports. To
be precise, let $\mathrm{supp}(\mu_{n})=\\{X_{1},\dots,X_{n}\\}$ and
$\mathrm{supp}(\nu_{m})=\\{Y_{1},\dots,Y_{m}\\}$. Define
$\mathcal{C}=\\{Z_{1},Z_{2},\dots,Z_{m+n}\\}=\\{X_{1},\dots,X_{n},Y_{1},\dots,Y_{m}\\}$.
Then discrete analog of equation 2.3 is

Formulation 1 (discrete):

$\displaystyle\begin{aligned} &\min_{\Pi,\mathbf{s}}&&\langle
C_{aug},\Pi\rangle+\lambda\left[\|s_{1}\|_{1}+\|t_{1}\|_{1}\right]\\\ &{\sf
subject\ to}&&\Pi 1_{m+n}=\begin{bmatrix}\mu_{n}+s_{1}\\\
t_{1}\end{bmatrix},\quad\Pi^{\top}1_{m+n}=\begin{bmatrix}0\\\
\nu_{m}\end{bmatrix}\\\ &&&\Pi\succeq 0,\quad
1_{m+n}^{\top}\mathbf{s}=0,\end{aligned}$ (2.6)

where $C_{aug}\in\mathbb{R}^{(m+n)\times(m+n)}$ is the augmented cost function
$C_{aug,i,j}=c(Z_{i},Z_{j})$ ($c$ is the ground cost, e.g., squared Euclidean
distance), $\mathbf{s}=(s_{1},t_{1})$ and $1_{r}$ is the vector all ones in
${\mathbb{R}}^{r}$. The TV norm got replaced with its discrete analog, the
$L_{1}$ norm. Similarly to its continuous counterpart, the optimization
problem is harder than typical OT due to additional constraint optimization
variable $\mathbf{s}$ and increased cost matrix size.

The discrete analog of equation 2.5 is straightforward:

Formulation 2 (discrete):

$\displaystyle\begin{aligned} &\min_{\Pi\in{\mathbb{R}}^{n\times m}}&&\langle
C_{\lambda},\Pi\rangle\\\ &{\sf subject\ to}&&\Pi
1_{n}=\mu_{n},\quad\Pi^{\top}1_{m}=\nu_{m},\quad\Pi\succeq 0,\end{aligned}$
(2.7)

where $C_{\lambda,i,j}=\min\\{c(X_{i},Y_{j}),2\lambda\\}$. As in the
continuous case, it is easy to adapt modern (regularized) OT solvers without
any computational overhead, and formulations of equation 2.6 and equation 2.7
are equivalent. It is also possible to recover $\mathbf{s}$ of equation 2.6
from the solution of equation 2.7 to perform outlier detection.

#### Two-sided formulation.

So far we have assumed that one of the input distributions does not have
outliers, which is the setting of MKE, where the clean distribution
corresponds to the model we learn. In some applications, both distributions
may be corrupted. To address this case, we provide an _equivalent_ two-sided
formulation, analogous to UOT with TV norm:

Formulation 3 (two-sided):

$\begin{split}\begin{aligned} \min_{\Pi,\mathbf{s}_{1},\mathbf{s}_{2}}&\ \ \
\langle
C_{aug},\Pi\rangle+\lambda\left[\|s_{1}\|_{1}+\|t_{1}\|_{1}+\|s_{2}\|_{1}+\|t_{2}\|_{1}\right]\\\
{\sf subject\ to}&\hskip 50.00008pt\Pi 1_{m+n}=\begin{bmatrix}\mu_{n}+s_{1}\\\
t_{1}\end{bmatrix}\\\ &\hskip
50.00008pt\Pi^{\top}1_{m+n}=\begin{bmatrix}s_{2}\\\
\nu_{m}+t_{2}\end{bmatrix}\\\ &\Pi\succeq 0,\quad
1_{m+n}^{\top}\mathbf{s}_{1}=0,\quad
1_{m+n}^{\top}\mathbf{s}_{2}=0\end{aligned}\end{split}$ (2.8)

where $\mathbf{s}_{1}=(s_{1}^{\top},t_{1}^{\top})^{\top}$ and
$\mathbf{s}_{2}=(s_{2}^{\top},t_{2}^{\top})^{\top}$.

##  3 Equivalence of the ROBOT formulations

In this section, we present our main theorem, which demonstrates the
equivalence between two formulations of the robust optimal transport:

######  Theorem 3.1.

For any two measures $\mu$ and $\nu$, ROBOT$(\mu,\nu)$ has same value for both
the formulations, i.e., Formulation 1 is equivalent to Formulation 2 for the
discrete case. Additionally, if the (non-truncated) cost function
$c(\cdot,\cdot)$ is a metric, then the equivalence of the two formulations
also holds for the continuous case. Moreover, we can recover optimal coupling
of one formulation from the other.

Below we sketch the proof of this theorem and highlight some important
techniques used in the proof. We focus on the discrete case as it is more
intuitive and has concrete practical implications in our experiments. A
complete proof can be found in Appendix A. Please also see Appendix A.2 for
the proof of equivalence between Formulations 1, 2 and 3 in the discrete case.

###  3.1 Proof sketch

In the remainder of this section, we consider the discrete case, i.e.,
equation 2.6 for Formulation 1 (F1) and equation 2.7 for Formulation 2 (F2).
Suppose $\Pi^{*}_{2}$ is an optimal solution of F2. Then we construct a
feasible solution $\Pi^{*}_{1},\mathbf{s}^{*}_{1}=(s^{*}_{1},t^{*}_{1})$ of F1
based on $\Pi^{*}_{2}$ with the same value of the objective function as F2 and
claim that $(\Pi^{*}_{1},\mathbf{s}^{*}_{1})$ is an optimal solution. We prove
the claim by contradiction: if $(\Pi^{*}_{1},\mathbf{s}^{*}_{1})$ is not
optimal, then there exists another pair
$(\tilde{\Pi}_{1},\tilde{\mathbf{s}}_{1})$ which is optimal for F1 with
strictly less objective value. We then construct another feasible solution
$\Pi^{*}_{2,new}$ of Formulation 2 which has the same objective value as of
$(\tilde{\Pi}_{1},\tilde{\mathbf{s}}_{1})$ for F1. This implies
$\Pi^{*}_{2,new}$ has strictly less objective value for F2 than $\Pi^{*}_{2}$,
which is a contradiction.

The two main steps of this proof are (1) constructing a feasible solution of
F1 starting from a feasible solution of F2 and (2) showing that the solution
constructed is indeed optimal for F1. Hence step (1) gives a recipe to
construct an optimal solution of F1 starting from an optimal solution of F2.
We elaborate the first point in the next subsection, which has practical
implications for outlier detection. The other point is more technical;
interested readers may go through the proof in Appendix A.1.

Algorithm 1 Generating optimal solution of F1 from F2

1: Start with $\Pi^{*}_{2}\in\mathbb{R}^{n\times m}$, an optimal solution of
Formulation 2.

2: Create an augmented matrix $\Pi\in\mathbb{R}^{m+n\times m+n}$ with all $0$.
Divide $\Pi$ into four blocks:
$\Pi=\begin{bmatrix}\underbrace{\Pi_{11}}_{n\times
n}&\underbrace{\Pi_{12}}_{n\times m}\\\ \underbrace{\Pi_{21}}_{m\times
n}&\underbrace{\Pi_{22}}_{m\times m}\end{bmatrix}$

3: Set $\Pi_{12}\leftarrow\Pi^{*}_{2}$ and collect all the indices
$\mathcal{I}=\\{(i,j):C_{i,j}>2\lambda\\}$.

4: Set $\Pi_{12}(i,j)\leftarrow 0$ for $(i,j)\in\mathcal{I}$.

5: Set
$\Pi_{22}(j,j)\leftarrow\sum_{i=1}^{n}\Pi^{*}_{2}(i,j)\mathds{1}_{(i,j)\in\mathcal{I}}$
for all $1\leq j\leq m$ and set $\Pi^{*}_{1}\leftarrow\Pi$.

6: Set
$s^{*}_{1}(i)=-\sum_{j=1}^{m}\Pi^{*}_{2}(i,j)\mathds{1}_{(i,j)\in\mathcal{I}}$
for all $1\leq i\leq n$.

7: Set $t_{1}^{*}(j)=\Pi_{22}(j,j)$ for all $1\leq j\leq m$.

8: return $\Pi^{*}_{1},s^{*}_{1},t^{*}_{1}$.

###  3.2 Going from Formulation 2 to Formulation 1

Let $\Pi^{*}_{2}$ (respectively $\Pi^{*}_{1}$) be an optimal solution of F2
(respectively F1). Recall that $\Pi_{1}^{*}$ has dimension $(m+n)\times(m+n)$.
From the column sum constraint in F1, we need to take the first $n$ columns of
$\Pi^{*}_{1}$ to be exactly $0$, whereas the last $m$ columns must sum up to
$\nu_{m}$. For any matrix $A$, we denote by $A[(a:b)\times(c:d)]$ the
submatrix consisting of rows from $a$ to $b$ and columns from $c$ to $d$. Our
main idea is to put a modified version of $\Pi^{*}_{2}$ in
$\Pi^{*}_{1}[(1:n)\times(n+1:m+n)]$ and make
$\Pi^{*}_{1}[(n+1:m+n)\times(n+1:m+n)]$ diagonal. First we describe how to
modify $\Pi^{*}_{2}$. Observe that, if for some $(i,j)$ $C_{i,j}>2\lambda$, we
expect $X_{i}\in\mathrm{supp}(\mu_{n})$ to be an outlier resulting in high
transportation cost, which is why we truncate the cost in F2. Therefore, to
get an optimal solution of F1, we make the corresponding value of optimal plan
$0$ and dump the mass into the corresponding slack variable $t_{1}^{*}$ in the
diagonal of the bottom right submatrix. This changes the row sum, which is
taken care of by $s_{1}^{*}$. But, as we are not moving this mass outside the
corresponding column, the column sum of $\Pi^{*}_{1}[(1:(m+n)):((n+1):(m+n))]$
remains same as column sum of $\Pi^{*}_{2}$, which is $\nu_{n}$. We summarize
this procedure in Algorithm 1.

Figure 1: Constructing optimal solution of Formulation 1 from optimal solution
of Formulation 2.

#### Example.

In Figure 1, we provide an example to visualize the construction. On the left,
we have $\Pi^{*}_{2}$, an optimal solution of Formulation 2. The blue
triangles denote the positions where the corresponding cost value is $\leq
2\lambda$, and light-green squares denote the positions where the
corresponding value of the cost matrix is $>2\lambda$. To construct an optimal
solution $\Pi^{*}_{1}$ of Formulation 1 from this $\Pi^{*}_{2}$, we first
create an augmented matrix of size $6\times 6$. We keep all the entries of the
left $6\times 3$ sub-matrix as $0$ (in this picture blank elements indicate
$0$). On the right submatrix, we put $\Pi^{*}_{2}$ into the top-right block,
but remove the masses from light-green squares, i.e. where cost value is
$>2\lambda$, and put it in the diagonal entries of the bottom right block as
shown in Figure 1. This mass contributes to the slack variables $s_{1}$ and
$t_{1}$, and this augmented matrix along with $s_{1},t_{1}$ give us an optimal
solution of Formulation 1.

###  3.3 Outlier detection with ROBOT

Our construction algorithm has practical consequences for outlier detection.
Suppose we have two datasets, a clean dataset $\nu_{m}$ (i.e., has no
outliers) and an outlier-contaminated dataset $\mu_{n}$. We can detect the
outliers in $\mu_{n}$ without directly solving costly Formulation 1 by
following Algorithm 2. In this algorithm, $\lambda$ is a regularization
parameter that can be chosen via cross-validation or heuristically (see
Section 4.2 for an example). In Section 4.2, we use this algorithm to perform
outlier detection on image data.

#### Outlier detection with entropic regularization.

Algorithm 1 allows us to recover solution of Formulation 1, which ultimately
is used for outlier detection in Algorithm 2, by solving the simpler truncated
cost Formulation 2 problem in equation 2.7. Similarly to the regular OT, it
can be solved exactly with a linear program solver—Pele & Werman (2009)
propose a faster exact solution based on min-cost-flow solvers benefiting from
the cost truncation—or approximately using entropic regularization techniques,
e.g. Sinkhorn algorithm (Cuturi, 2013). In the following lemma, we show that
Algorithm 1 recovers a meaningful approximate solution of Formulation 1 from
an approximate solution of Formulation 2 obtained with entropy-regularized OT
solvers.

######  Lemma 3.2.

Let $\Pi_{2,\alpha}^{*}$ be a solution of the entropy regularized version of
equation 2.7:

$\displaystyle\begin{aligned}
&\operatorname*{arg\,min}_{\Pi\in{\mathbb{R}}^{n\times m}}&&\langle
C_{\lambda},\Pi\rangle+\alpha H(\Pi)\\\ &{\sf subject\ to}&&\Pi
1_{n}=\mu_{n},\quad\Pi^{\top}1_{m}=\nu_{m},\quad\Pi\succeq 0,\end{aligned}$
(3.1)

and let $\left(\Pi_{1,\alpha}^{*},\mathbf{s}^{*}_{1,\alpha}\right)$ be the
corresponding approximate solution of Formulation 1 recovered from
$\Pi^{*}_{2,\alpha}$ by Algorithm 1. Then

$\|\Pi^{*}_{1,\alpha}-\Pi^{*}_{1}\|_{F}+\|\mathbf{s}^{*}_{1,\alpha}-\mathbf{s}^{*}_{1}\|_{2}\rightarrow
0$

as $\alpha\to 0$, where $\left(\Pi^{*}_{1},\mathbf{s}^{*}_{1}\right)$ is the
exact solution of Formulation 1 in equation 2.6.

When the solution of equation 2.7 is non-unique, then the solution of equation
3.1 converges to the solution of equation 2.7 with maximum entropy (see
Proposition 4.1 of Peyré & Cuturi (2018)) and consequently recovers the
corresponding maximum entropy solution
$\left(\Pi_{1}^{*},\mathbf{s}^{*}_{1}\right)$ of equation 2.6 by Algorithm 1
in the limit. Please find the proof in Appendix C.

Algorithm 2 Outlier detection in contaminated data

1: Start with $\mu_{n}$ (contaminted data) and $\nu_{m}$ (clean data).

2: Solve Formulation 2 and obtain $\Pi^{*}_{2}$ using a suitable value of
$\lambda$.

3: Use Algorithm 1 to obtain $\Pi^{*}_{1},s^{*}_{1},t^{*}_{1}$ from
$\Pi^{*}_{2}$.

4: Find $\mathcal{I}$, the set of all the indices where $\mu_{n}+s_{1}^{*}=0$.

5: Return $\mathcal{I}$ as the indices of outliers in $\mu_{n}$.

##  4 Empirical studies

To evaluate effectiveness of ROBOT, we consider the task of robust mean
estimation under the Huber contamination model. The data is generated from
$(1-{\epsilon})\mathcal{N}(\eta_{0},I_{d})+{\epsilon}\mathcal{N}(\eta_{1},I_{d})$
and the goal is to estimate $\eta_{0}$. Prior work has advocated for using
$f$-divergence GANs (Chao et al., 2018; Wu et al., 2020) for this problem and
pointed out inefficiencies of Wasserstein GAN in the presence of outliers. We
show that our robust OT formulation allows to estimate the uncontaminated mean
$\eta_{0}$ comparably or better than a variety of $f$-divergence GANs. We also
use this simulated setup to study sensitivity to the cost truncation
hyperparameter $\lambda$.

In our second experiment, we present a new application of optimal transport
enabled by ROBOT. Suppose we have collected a curated dataset $\nu_{m}$ (i.e.,
we know that it has no outliers)—such data collection is expensive, and we
want to benefit from it to automate subsequent data collection. Let $\mu_{n}$
be a second dataset collected “in the wild,” i.e., it may or may not have
outliers. We demonstrate how ROBOT can be used to identify outliers in
$\mu_{n}$ using the curated dataset $\nu_{m}$.

Table 1: Robust mean estimation with GANs using different distribution divergences. True mean is $\eta_{0}=\mathbf{0}_{5}$; sample size $n=1000$; contamination proportion ${\epsilon}=0.2$. We report results over 30 experiment restarts. Contamination | JS Loss | SH Loss | RKL Loss | ROBOT | UOT
---|---|---|---|---|---
$\mathcal{N}(0.1\cdot\mathbf{1_{5}},I_{5})$ |  0.09 $\pm$ 0.03 | 0.11 $\pm$ 0.03 | 0.115 $\pm$ 0.03 | 0.1 $\pm$ 0.03 | 0.1 $\pm$ 0.04
$\mathcal{N}(0.5\cdot\mathbf{1_{5}},I_{5})$ | 0.23 $\pm$ 0.04 | 0.24 $\pm$ 0.05 | 0.24 $\pm$ 0.05 |  0.117 $\pm$ 0.03 | 0.2 $\pm$ 0.04
$\mathcal{N}(1\cdot\mathbf{1_{5}},I_{5})$ | 0.43 $\pm$ 0.05 | 0.43 $\pm$ 0.06 | 0.43 $\pm$ 0.06 | 0.261 $\pm$ 0.06 |  0.25 $\pm$ 0.05
$\mathcal{N}(2\cdot\mathbf{1_{5}},I_{5})$ | 0.67 $\pm$ 0.07 | 0.67 $\pm$ 0.08 | 0.67 $\pm$ 0.08 | 0.106 $\pm$ 0.03 |  0.1 $\pm$ 0.03

(a) Varying proportion of contamination (b) Varying outlier distribution mean

Figure 2: Empirical study of the cost truncation hyperparameter $\lambda$
sensitivity.

###  4.1 Robust mean estimation

Following Wu et al. (2020), we consider a simple generator of the form
$g_{\theta}(x)=x+\theta$, $x\sim\mathcal{N}(0,I_{d})$, $d$ is the data
dimension. The basic idea of robust mean estimation with GANs is to minimize
various distributional divergences between samples from $g_{\theta}$ and
observed data simulated from
$(1-{\epsilon})\mathcal{N}(\eta_{0},I_{d})+{\epsilon}\mathcal{N}(\eta_{1},I_{d})$.
The goal is to estimate $\eta_{0}$ with $\theta$.

To efficiently implement ROBOT GAN, we use a standard min-max optimization
approach: solve the inner max (ROBOT) and use gradient descent for the outer
min parameter. To solve ROBOT, it is straightforward to adopt any of the prior
stochastic regularized OT solvers: the only modification is the truncation of
the cost entries as in equation 2.7. We use the stochastic algorithm for semi-
discrete regularized OT (Genevay et al., 2016, Algorithm 2). We summarize
ROBOT GAN in Algorithm 3. Line 5 - Line 11 perform the inner optimization
where we solve the entropy regularized OT dual with truncated cost and Line 12
- Line 14 perform gradient update of $\theta$.

For the $f$-divergence GANs (Nowozin et al., 2016), we use the code of Wu et
al. (2020) for GANs with Jensen-Shannon (JS) loss, squared Hellinger (SH)
loss, and Reverse Kullback-Leibler (RKL) loss. For the exact expressions of
these divergences, see Table 1 of Wu et al. (2020). We report estimation error
measured by the Euclidean distance between true uncontaminated mean $\eta_{0}$
and estimated mean $\theta$ for various contamination distributions in Table
1. ROBOT GAN performs well across all considered contamination distributions.
As the difference between true mean $\eta_{0}$ and contamination mean
$\eta_{1}$ increases, the estimation error of all methods tends to increase.
However, when it becomes easier to distinguish outliers from clean samples,
i.e., $\eta_{1}=2\cdot\mathbf{1_{5}}$, performance of ROBOT noticeably
improves. We present an analogous study with data simulated from a mixture of
Cauchy distributions in Appendix F111All codes are available at
https://github.com/debarghya-mukherjee/Robust-Optimal-Transport..

Algorithm 3 ROBOT GAN

1: Input:  robustness regularizion $\lambda$, entropic regularization
$\alpha$, data distribution $\mu_{n}\in\Delta^{n-1}$,
$supp(\mu_{n})=\mathcal{X}=[X_{1},\dots,X_{n}]$, steps sizes $\tau$ and
$\gamma$

2: Initialize:  Initialize $\theta=\theta_{init}$, set number of iterations
$M$ and $L$, $i=0$, $\mathbf{v}=\tilde{\mathbf{v}}=\mathbf{0}$.

3: for $j=1,\dots,M$ do

4: Generate $\tilde{z}\sim\mathcal{N}(0,I_{d})$ and set $z=\tilde{z}+\theta$.

5: Set the cost vector $\mathbf{c}\in\mathbb{R}^{n}$ as
$\mathbf{c}(k)=\min\\{c(X_{k},z),2\lambda\\}$ for $k=1,\dots,n$.

6: for $i=1,\dots,L$ do

7: Set $\mathbf{h}\leftarrow\frac{\tilde{\mathbf{v}}-\mathbf{c}}{\alpha}$ and
do the normalized exponential transformation
$\mathbf{u}\leftarrow\frac{e^{\mathbf{h}}}{\langle\mathbf{1},e^{\mathbf{h}}\rangle}$.

8: Calculate the gradient
$\nabla\tilde{\mathbf{v}}\leftarrow\mu_{n}-\mathbf{u}$.

9: Update
$\tilde{\mathbf{v}}\leftarrow\tilde{\mathbf{v}}+\gamma\nabla\tilde{\mathbf{v}}$
and $\mathbf{v}\leftarrow(1/(j+i))\tilde{\mathbf{v}}+(j+i-1/(j+i))\mathbf{v}$.

10: end for

11: Do the same transformation of $\mathbf{v}$ as in Step 7, i.e. set
$\mathbf{h}\leftarrow\frac{\mathbf{v}-\mathbf{c}}{\alpha}$ and set
$\Pi\leftarrow\frac{e^{\mathbf{h}}}{\langle\mathbf{1},e^{\mathbf{h}}\rangle}$.

12: Set $\Pi(k)=0$ for $k$ such that $c(X_{k},z)>2\lambda$ for $k=1,\dots,n.$

13: Calculate gradient with respect to $\theta$ as
$\nabla\theta=2\left[z\sum_{k}\Pi(k)-\mathcal{X}^{\top}\Pi\right]$

14: Update $\theta\leftarrow\theta-\tau\nabla\theta$.

15: end for

16: Ouput: $\theta$

We also compared to the Sinkhorn-based UOT algorithm (Chizat et al., 2018)
available in the Python Optimal Transport (POT) library (Flamary & Courty,
2017); to obtain a UOT GAN, we modified steps 5-12 of Algorithm 3 for
computing $\Pi$. Unsurprisingly, both ROBOT and UOT perform similarly: recall
equivalence to Formulation 3, which is similar to UOT with TV norm. The key
insight of our work is the equivalence to classical OT with truncated cost,
that greatly simplifies optimization and allows to use existing stochastic OT
algorithms. In this experiment, the sample size $n=1000$ is sufficiently small
for the Sinkhorn-based UOT POT implementation to be effective, but it breaks
in the experiment we present in Section 4.2. We also tried the code of Balaji
et al. (2020) based on CVXPY (Diamond & Boyd, 2016), but it is too slow even
for the $n=1000$ sample size. In Subsection 4.3 we present a comparison to
Balaji et al. (2020) on a smaller sample size.

#### Hyperparameter sensitivity study.

In the previous experiment, we set $\lambda=0.5$. Now we demonstrate
empirically that there is a broad range of $\lambda$ values performing well.
In Figure 2(a), we study sensitivity of $\lambda$ under various contamination
proportions ${\epsilon}$ holding $\eta_{0}=\mathbf{1}_{5}$ and
$\eta_{1}=5\cdot\mathbf{1}_{5}$ fixed. Horizontal lines correspond to
$\lambda=\infty$, i.e., vanilla OT. The key observations are: there is a wide
range of $\lambda$ efficient at all contamination proportions (note the
$\log_{10}$ $x$-axis scale), and ROBOT is always at least as good as vanilla
OT (even when there is no contamination ${\epsilon}=0$). In Figure 2(b), we
present a similar study varying the mean of the contamination distribution and
holding ${\epsilon}=0.2$ fixed. We see that as the contamination distribution
gets closer to the true distribution, it becomes harder to pick a good
$\lambda$, but the performance is always at least as good as the vanilla OT
(horizontal lines).

###  4.2 Outlier detection for data collection

Our robust OT formulation equation 2.6 enables outlier identification. Let
$\nu_{m}$ be a clean dataset and $\mu_{n}$ potentially contaminated with
outliers. Recall that ROBOT allows modification of one of the input
distributions to eliminate potential outliers. We can identify outliers in
$\mu_{n}$ as follows: if $\mu_{n}(i)+s^{*}_{1}(i)=0$, then $X_{i}$, the $i$th
point in $\mu_{n}$, is an outlier. Instead of directly solving equation 2.6,
which may be inefficient, we use our equivalence results and solve an easier
optimization problem equation 2.7, followed by recovering $\mathbf{s}$ to find
outliers via Algorithm 2. When using entropy-regularized approximate solutions
to detect outliers with Algorithm 2, in step 4, $\mu_{n}+s^{*}_{1}$ will not
be exactly 0 for the outliers, so a small threshold should be used instead. We
modify step 4 to “Find $\mathcal{I}$, the set of all the indices where
$\mu_{n}+s_{1}^{*}<1/n^{2}$” when using entropy regularization.

To test our outlier-identification methodology we follow the experimental
setup of Tagasovska & Lopez-Paz (2019). Specifically, let $\nu_{m}$ be a clean
dataset consisting of 10k MNIST digits from 0 to 4 and $\mu_{n}$ be a dataset
collected “in the wild” consisting of (different) 8k MNIST digits from 0 to 4
and 2k outlier MNIST images of digits from 5 to 9. We compute
ROBOT$(\mu_{n},\nu_{m})$ to identify outlier digit images in $\mu_{n}$. For
each point in $\mu_{n}$ we obtain a prediction, outlier or clean, which allows
us to evaluate accuracy. Tagasovska & Lopez-Paz (2019) use last-layer features
of a neural network pre-trained on the clean data $\nu_{m}$—it is
straightforward to combine ROBOT and other baselines we consider with any
feature extractor, but in this experiment we simply use the raw data.

We compare to the Orthonormal Certificates (OC) method of Tagasovska & Lopez-
Paz (2019) and to a variety of standard outlier detection algorithms available
in Scikit-learn (Pedregosa et al., 2011): one class SVM (Schölkopf et al.,
1999), local outlier factor (Breunig et al., 2000), isolation forest (Liu et
al., 2008) and elliptical envelope (Rousseeuw & Driessen, 1999). All baselines
except one class SVM and local outlier factor use clean data for training as
does our method.

Table 2: Outlier detection on MNIST. Methods | Accuracy
---|---
One class SVM | $0.496\pm 0.003$
Local outlier factor | $0.791\pm 0.001$
Isolation forest | $0.636\pm 0.010$
Elliptical envelope | $0.739\pm 0.002$
Orthonormal Certificates | $0.819\pm 0.008$
ROBOT | $0.859\pm 0.002$
ROBOT-Sinkhorn | $\textbf{0.897}\pm 0.004$

$\overbrace{\rule{95.00014pt}{0.0pt}}^{\text{Inliers detected as outliers}}$ | $\overbrace{\rule{95.00014pt}{0.0pt}}^{\text{Outliers detected as inliers}}$
---|---
digit | count | digit | count
0 | 10 | 5 | 191
1 | 1 | 6 | 94
2 | 53 | 7 | 175
3 | 32 | 8 | 181
4 | 22 | 9 | 309

[table]A table beside a figure

Figure 3: Insights into ROBOT-Sinkhorn performance: the top left figure is a
collection of random inlier digits miss-classified as outliers; the top right
picture represents random outlier digits miss-classified as inliers and the
table illustrates frequency distribution of miss-classified images. The
majority of the errors are on digit 9, possibly due to its similarity to 4.

Results of 30 experiment repetitions are summarized in Table 2. ROBOT, i.e.
Algorithm 2 where equation 2.7 is solved exactly with a linear program solver,
and ROBOT-Sinkhorn, i.e. Algorithm 2 where equation 2.7 is solved
approximately using Sinkhorn (Cuturi, 2013), produce the best results. In this
experiment ROBOT-Sinkhorn outperforms ROBOT, but it is not necessarily to be
expected in general. We conclude that our method is effective in assisting
data collection once an initial set of clean data has been acquired and is
compatible with entropy-regularized OT solvers ensuring scalability.

Elliptical envelope assumes that clean data is Gaussian, while one class SVM,
isolation forest and local outlier factor correspondingly attempt to fit SVM,
random forest and $k$-nearest neighbors classifiers to distinguish clean and
outlier samples. All these baselines work best when the clean data is
unimodal, which is not the case for the MNIST 0 to 4 digits considered in our
experiment. The orthonormal certificates method is rooted in PCA and assumes
that clean and outlier data live in different subspaces. This assumption is
reasonable for the MNIST data, but it might not hold broadly in practice. We
believe that empirical success of our method is due to its optimal transport
nature. OT is a geometry-sensitive metric on distributions that can
distinguish multi-modal distributions and distributions supported on the same
subspace. We provide additional insights into the ROBOT performance in Figure
3. A theoretical investigation of ROBOT outlier-detection guarantees is an
interesting future work direction.

#### Hyperparameter selection.

To select the cost truncation hyperparameter $\lambda$, we propose the
following heuristic: since we know that $\nu_{m}$ is clean, we can subsample
two datasets from it, compute vanilla OT to obtain transportation plan $\Pi$
and set $\lambda$ to be half the maximum distance between matched elements,
i.e. $2\lambda=\max_{i,j}\\{C_{ij}:\,\Pi_{ij}>0\\}$, where $C$ is the cost
matrix for the two subsampled datasets. This procedure is essentially
estimating maximum distance between matched clean samples. To avoid
subsampling noise we use 99th percentile instead of the maximum. Our
experiments also revealed that increasing $\lambda$ increases the set of
outliers progressively, i.e. if a sample is detected as outlier for some value
of $\lambda$, then it will be detected as outlier for all higher values of
$\lambda$. A rigorous theoretical analysis for this observation is a potential
future work.

The Orthonormal Certificates method (Tagasovska & Lopez-Paz, 2019) also
requires setting a threshold. It computes the null-space of the clean train
data and uses the norm of the projection into that space as a score to
distinguish outliers. Similarly to ROBOT, and as the authors do in their code,
we use 99th percentile of those scores computed on the clean train data as the
threshold. For other baselines, we use the default hyperparameters.

###  4.3 Comparison with Balaji et al. (2020)

We conduct additional experiments comparing to the recent robust optimal
transport method of Balaji et al. (2020). Their method relies on CVXPY and
does not scale to the sample sizes considered in our previous experiments. We
compare on smaller data sizes.

Robust mean estimation. We set $n=200$ samples with contamination distribution
mean equal to 2 and true mean equal to $0$ (same configuration as in the last
row in Table 1). The runtime and the estimation error is reported in Table 3.

Table 3: Robust mean estimation for $n=200$. Method | Estimation error | Run-time
---|---|---
ROBOT-Sinkhorn | $0.196\pm 0.1$ | $35s\pm 0.6s$
Balaji et al. (2020) | $0.212\pm 0.074$ | $7745s\pm 1670s$

Outlier detection experiment. In this experiment, we consider the size of the
entire dataset (inliers + outliers) to be $n=1000$ (with $800$ inliers and
$200$ outliers). The results (accruacy and run-time) are provided in Table 4.

Table 4: Outlier detection for $n=1000$. Method | Estimation error | Run-time
---|---|---
ROBOT-Sinkhorn | $0.86\pm 0.015$ | $6\pm 2s$
Balaji et al. (2020) | $0.8\pm 0.0005$ | $3343\pm 960s$

The method of Balaji et al. (2020) is significantly slower as expected.
Comparing the performance, we think that ROBOT performs better because it is
based on TV norm, while the method of Balaji et al. (2020) uses a chi-squared
constraints on the marginal perturbations. A TV constraint on the marginal
perturbations is more closely related to the ${\epsilon}$-contamination model
for outlier detection, suggesting that a TV-based constraint/regularizer could
be a better choice for the outlier detection applications. We also note that
in the outlier detection experiment, using chi-square divergence results in a
non-sparse solution and requires tuning a threshold parameter to perform
outlier detection, in addition to the chi-square distance hyperparameter. We
tuned those parameters, but were not able to achieve significant performance
improvements for the method of Balaji et al. (2020).

##  5 Summary and discussion

We propose and study ROBOT, a robust formulation of optimal transport.
Although the problem is seemingly asymmetric and challenging to optimize,
there is an equivalent formulation based on cost truncation that is symmetric
and compatible with modern stochastic optimization methods for OT.

ROBOT closely resembles unbalanced optimal transport (UOT). In our
formulation, we added a TV regularizer to the vanilla optimal transport
problem. This is motivated by the ${\epsilon}$-contamination model. In UOT,
the TV regularizer is typically replaced with a KL divergence. Other choices
of the regularizer may lead to new properties and applications. Studying
equivalent, simpler formulations of UOT with different divergences may be a
fruitful future work direction.

From the practical perspective, in our experiments we observed no degradation
of ROBOT GAN in comparison to OT GAN, even when there were no outliers. It is
possible that replacing OT with ROBOT may be beneficial for various machine
learning applications of OT. Data encountered in practice may not be
explicitly contaminated with outliers, but it often has errors and other
deficiencies, suggesting that a “no-harm” robustness is desirable.

## Acknowledgements

This paper is based upon work supported by the National Science Foundation
(NSF) under grants no. 1830247 and 1916271. J. Solomon acknowledges the
generous support of Army Research Office grants W911NF1710068 and
W911NF2010168, of Air Force Office of Scientific Research award
FA9550-19-1-031, of National Science Foundation grant IIS-1838071, from the
CSAIL Systems that Learn program, from the MIT–IBM Watson AI Laboratory, from
the Toyota–CSAIL Joint Research Center, from a gift from Adobe Systems, and
from the Skoltech–MIT Next Generation Program. We also thank anonymous
reviewers, whose comments were extremely helpful for further improvement of
our paper.

##  Appendix A Proof of Theorem 3.1

In the proofs, $a\wedge b$ denotes $\min\\{a,b\\}$ for any
$a,b\in{\mathbb{R}}$.

###  A.1 Proof of discrete version

###### Proof.

Define a matrix $\Pi$ as:

$\Pi(i,j)=\begin{cases}0,&\text{if }C(i,j)>2\lambda\\\
\Pi^{*}_{2}(i,j),&\text{otherwise}\end{cases}$

Also define $s\in\mathbb{R}^{n}$ and $t\in\mathbb{R}^{m}$ as:

$s^{*}_{1}(i)=-\sum_{j=1}^{m}\Pi^{*}_{2}(i,j)\mathds{1}_{C(i,j)>2\lambda}$

and similarly define:

$t^{*}_{1}(j)=\sum_{i=1}^{n}\Pi^{*}_{2}(i,j)\mathds{1}_{C(i,j)>2\lambda}$

These vectors corresponds to the row sums and the column sums of the elements
of the optimal transport plan of Formulation 2, where the cost function
exceeds $2\lambda$. Note that, these co-ordinates of the optimal transport
plan corresponding to those co-ordinates of cost matrix, where the cost is
greater than $2\lambda$ and contribute to the objective value via their sum
only, hence any different arrangement of these transition probabilities with
same sum gives the same objective value.

Now based on this $\Pi$ obtained we construct a feasible solution of
Formulation 1 following Algorithm 1:

$\Pi^{*}_{1}=\begin{bmatrix}\mathbf{0}&\Pi\\\ \mathbf{0}&{\sf
diag}(t^{*}_{1})\end{bmatrix}$

The row sums of $\Pi^{*}_{1}$ is:

$\Pi^{*}_{1}\mathbf{1}=\begin{bmatrix}\mu_{n}+s^{*}_{1}\\\
t^{*}_{1}\end{bmatrix}$

and it is immediate from the construction that the column sums of
$\Pi^{*}_{1}$ is $\nu_{m}$. Also as:

$\sum_{i=1}^{n}s^{*}_{1}(i)=\sum_{j=1}^{m}t^{*}_{1}(j)=\sum_{(i,j):C_{i,j}>2\lambda}\Pi^{*}_{2}(i,j)$

and $s^{*}_{1}\preceq 0,t^{*}_{1}\succeq 0$, we have:

$\mathbf{1}^{\top}(\mu_{n}+s^{*}_{1}+t^{*}_{1})=\mathbf{1}^{\top}p=1\,.$

Therefore, we have $(\Pi^{*}_{1},s^{*}_{1},t^{*}_{1})$ is a feasible solution
of Formulation 1. Now suppose this is not an optimal solution. Pick an optimal
solution $\tilde{\Pi},\tilde{s},\tilde{t}$ of Formulation 1 so that:

$\langle
C_{aug},\tilde{\Pi}\rangle+\lambda\left[\|\tilde{s}\|_{1}+\|\tilde{t}\|_{1}\right]<\langle
C_{aug},\Pi^{*}_{1}\rangle+\lambda\left[\|s^{*}_{1}\|_{1}+\|t^{*}_{1}\|_{1}\right]$

The following two lemmas provide some structural properties of any optimal
solution of Formulation 1:

######  Lemma A.1.

Suppose $\Pi^{*}_{1},s^{*}_{1},t^{*}_{1}$ are optimal solution for Formulation
1. Divide $\Pi^{*}_{1}$ into four parts corresponding to augmentation as in
algorithm 1:

$\Pi^{*}_{1}=\begin{bmatrix}\Pi^{*}_{1,11}&\Pi^{*}_{1,12}\\\
\Pi^{*}_{1,21}&\Pi^{*}_{1,22}\end{bmatrix}$

Then we have $\Pi^{*}_{1,11}=\Pi^{*}_{1,21}=\mathbf{0}$ and $\Pi^{*}_{1,22}$
is a diagonal matrix.

######  Lemma A.2.

If $\Pi^{*}_{1},s^{*}_{1},t^{*}_{1}$ is an optimal solution of Formulation 1
then:

  1. 1.

If $C_{i,j}>2\lambda$ then $\Pi^{*}_{1}(i,j)=0$.

  2. 2.

If $C_{i,j}<2\lambda$ for some $i$ and for all $1\leq j\leq n$, then
$s^{*}_{1}(i)=0$.

  3. 3.

If $C_{i,j}<2\lambda$ for some $j$ and for all $1\leq i\leq m$, then
$t^{*}_{1}(j)=0$.

  4. 4.

If $C_{i,j}<2\lambda$ then $s^{*}_{1}(i)t^{*}_{1}(j)=0$.

We provide the proofs in the next subsection. By Lemma A.1 we can assume
without loss of generality:

$\tilde{\Pi}=\begin{bmatrix}\mathbf{0}&\tilde{\Pi}_{12}\\\ \mathbf{0}&{\sf
diag}(\tilde{t})\end{bmatrix}$

Now based on $\left(\tilde{\Pi},\tilde{s},\tilde{t}\right)$ we create a
feasible solution namely $\Pi^{*}_{2,new}$ of Formulation 2 as follows: Define
the set of indices $\\{i_{1},\cdots,i_{k}\\}$ and $\\{j_{1},\dots,j_{l}\\}$
as:

$\tilde{s}_{i_{1}},\tilde{s}_{i_{2}},\dots,\tilde{s}_{i_{k}}>0\ \ \ \
\text{and}\ \ \
\tilde{t}_{j_{1}},\tilde{t}_{j_{2}},\dots,\tilde{t}_{j_{l}}>0\,.$

Then by part (4) of Lemma A.2 we have $C_{i_{\alpha},j_{\beta}}>2\lambda$ for
$\alpha\in\\{1,\dots,k\\}$ and $\beta\in\\{1,\dots,l\\}$. Also by part (2) of
Lemma A.2 the value of transport plan at these co-ordinates is 0. Now
distribute the mass of slack variables in these co-ordinates such that the
marginals of new transport plan becomes exactly $\mu_{n}$ and $\nu_{m}$. This
new transport plan is our $\Pi^{*}_{2,new}$. Recall that,
$\|\tilde{s}\|_{1}=\|\tilde{t}\|_{1}$. Hence, here the regularizer value
decreases by $2\lambda\|\tilde{s}\|_{1}$ and the cost value increased by
exactly $2\lambda\|\tilde{s}\|_{1}$ as we are truncating the cost. Hence we
have:

$\displaystyle\langle C_{\lambda},\Pi^{*}_{2,new}\rangle$
$\displaystyle=\langle
C_{aug},\tilde{\Pi}\rangle+\lambda\left[\|\tilde{s}\|_{1}+\|\tilde{t}\|_{1}\right]$
$\displaystyle<\langle
C_{aug},\Pi^{*}_{1}\rangle+\lambda\left[\|s^{*}_{1}\|_{1}+\|t^{*}_{1}\|_{1}\right]$
$\displaystyle=\langle C_{\lambda},\Pi^{*}_{2}\rangle$

which is contradiction as $\Pi^{*}_{2}$ is the optimal solution of Formulation
2. This completes the proof for the discrete part.

∎

###  A.2 Proof of equivalence for two sided formulation

Here we prove that our two sided formulation, i.e. Formulation 3 (equation
2.8) is equivalent to Formulation 1 (equation 2.6) for the discrete case.
Towards that end, we introduce another auxiliary formulation and show that
both Formulation 1 and Formulation 3 are equivalent to the following auxiliary
formulation of the problem.

Formulation 4:

$\begin{aligned} &\min\nolimits_{\Pi\in{\mathbb{R}}^{m\times
n},s_{1}\in{\mathbb{R}}^{m},s_{2}\in{\mathbb{R}}^{n}}&&\langle
C,\Pi\rangle+\lambda\left[\|s_{1}\|_{1}+\|s_{2}\|_{1}\right]\\\ &{\sf subject\
to}&&\Pi 1_{n}=p+s_{1}\\\ &&&\Pi^{T}1_{m}=q+s_{2}\\\ &&&\Pi\succeq
0\end{aligned}.$ (A.1)

First we show that Formulation 1 and Formulation 4 are equivalent in a sense
that they have the same optimal objective value.

######  Theorem A.3.

Suppose $C$ is a cost function such that $C(x,x)=0$. Then Formulation 1 and
Formulation 4 has same optimal objective value.

###### Proof.

Towards that end, we show that given one optimal variables of one formulation
we can get optimal variables of other formulation with the same objective
value. Before going into details we need the following lemma whose proof is
provided in Appendix B:

######  Lemma A.4.

Suppose $\Pi^{*}_{4},s^{*}_{4,1},s^{*}_{4,2}$ are the optimal variables of
Formulation 4. Then $s^{*}_{4,1}\preceq 0$ and $s^{*}_{4,2}\preceq 0$.

Now we prove that optimal value of Formulation 1 and Formulation 4 are same.
Let $(\Pi^{*}_{1},s^{*}_{1,1},t^{*}_{1,1})$ is an optimal solution of
Formulation 1. Then we claim that $(\Pi^{*}_{1},s^{*}_{1,1},t^{*}_{1,1})$ is
also an optimal solution of Formulation 4. Clearly it is feasible solution of
Formulation 4. Suppose it is not optimal, i.e. there exists another optimal
solution $(\tilde{\Pi}_{4},\tilde{s}_{4,1},\tilde{s}_{4,2})$ such that:

$\langle
C,\tilde{\Pi}_{4}\rangle+\lambda(\|\tilde{s}_{4,1}\|_{1}+\|\tilde{s}_{4,2}\|_{2})<\langle
C,\Pi^{*}_{1,12}\rangle+\lambda(\|s^{*}_{1,1}\|_{1}+\|t^{*}_{1,1}\|_{1})$

Now based on $(\tilde{\Pi}_{4},\tilde{s}_{4,1},\tilde{s}_{4,2})$ we construct
a feasible solution of Formulation 1 as follows:

$\tilde{\Pi}_{1}=\begin{bmatrix}\mathbf{0}&\tilde{\Pi}_{4}\\\ \mathbf{0}&-{\sf
diag}(\tilde{s}_{4,2})\end{bmatrix}$

Note that we proved in Lemma A.4 $\tilde{s}_{4,2}\preceq 0$, hence we have
$\tilde{\Pi}_{1}\succeq 0$. Now as the column sums of $\tilde{\Pi}_{4}$ is
$q+\tilde{s}_{4,2}$, we have column sums of $\tilde{\Pi}_{1}=[\mathbf{0}\ \
q^{\top}]^{\top}$ and the row sums are $[(p+\tilde{s}_{4,1})^{\top}\ \ \
\tilde{s}_{4,2}^{\top}]^{\top}$. Hence we take
$\tilde{s}_{1,1}=\tilde{s}_{4,1}$ and $\tilde{s}_{1,2}=\tilde{s}_{4,2}$. Then
it follows:

$\displaystyle\langle
C_{aug},\tilde{\Pi}_{1}\rangle+\lambda\left[\|\tilde{s}_{1,1}\|_{1}+\|\tilde{s}_{1,2}\|_{1}\right]$
$\displaystyle=\langle
C,\tilde{\Pi}_{4}\rangle+\lambda\left[\|\tilde{s}_{4,1}\|_{1}+\|\tilde{s}_{4,2}\|_{1}\right]$
$\displaystyle<\langle
C,\Pi^{*}_{1,12}\rangle+\lambda\left[\|s^{*}_{1,1}\|_{1}+\|t^{*}_{1,1}\|_{1}\right]$
$\displaystyle=\langle
C_{aug},\Pi^{*}_{1}\rangle+\lambda\left[\|s^{*}_{1,1}\|_{1}+\|t^{*}_{1,1}\|_{1}\right]$

This is contradiction as we assumed $(\Pi^{*}_{1},s^{*}_{1,1},t^{*}_{1,2})$ is
an optimal solution of Formulation 1. Therefore we conclude
$(\Pi^{*}_{1},s^{*}_{1,1},t^{*}_{1,1})$ is also an optimal solution of
Formulation 4 which further concludes Formulation 1 and Formulation 4 have
same optimal values. This completes the proof of the theorem. ∎

######  Theorem A.5.

The optimal objective value of Formulation 3 and Formulation 4 are same.

###### Proof.

Like in the proof of Theorem A.3 we also prove couple of lemmas.

######  Lemma A.6.

Any optimal transport plan $\Pi^{*}_{3}$ of Formulation 3 has the following
structure: If we write,

$\Pi^{*}_{3}=\begin{bmatrix}\Pi^{*}_{3,11}&\Pi^{*}_{3,12}\\\
\Pi^{*}_{3,21}&\Pi^{*}_{3,22}\end{bmatrix}$

then $\Pi^{*}_{3,11}$ and $\Pi^{*}_{3,22}$ are diagonal matrices and
$\Pi^{*}_{3,21}=\mathbf{0}$.

######  Lemma A.7.

If $s^{*}_{3,1},t^{*}_{3,1},s^{*}_{3,2},t^{*}_{3,2}$ are four optimal slack
variables in Formulation 3, then $s^{*}_{3,1},t^{*}_{3,1}\preceq 0$ and
$s^{*}_{3,2},t^{*}_{3,2}\succeq 0$.

###### Proof.

The line of argument is same as in proof of Lemma A.4. ∎

Next we establish equivalence. Suppose
$(\Pi^{*}_{3},s^{*}_{3,1},t^{*}_{3,1},s^{*}_{3,2},t^{*}_{3,2})$ are optimal
values of Formulation 3. We claim that
$(\Pi^{*}_{3,12},s^{*}_{3,1}-s^{*}_{3,2},t^{*}_{3,1}-t^{*}_{3,2})$ forms an
optimal solution of Formulation 4. The objective value will then also be same
as $s^{*}_{3,1}\preceq 0,s^{*}_{3,2}\succeq 0$ (Lemma A.7) implies
$\|s^{*}_{3,1}-s^{*}_{3,2}\|_{1}=\|s^{*}_{3,1}\|_{1}+\|s^{*}_{3,2}\|_{1}$ and
similarly $t^{*}_{3,1}\preceq 0,t^{*}_{3,2}\succeq 0$ implies
$\|t^{*}_{3,1}-t^{*}_{3,2}\|_{1}=\|t^{*}_{3,1}\|_{1}+\|t^{*}_{3,2}\|_{1}$.
Feasibility is immediate. Now for optimality, we again prove by contradiction.
Suppose they are not optimal. Then lets say
$\tilde{\Pi}_{4},\tilde{s}_{4,1},\tilde{s}_{4,2}$ are an optimal triplet of
Formulation 4. Now construct another feasible solution of Formulation 3 as
follows: Set
$\tilde{s}_{3,2}=\tilde{t}_{3,2}=0,\tilde{s}_{3,1}=\tilde{s}_{4,1}$ and
$\tilde{t}_{3,1}=\tilde{s}_{4,2}$. Set the matrix as:

$\tilde{\Pi}_{3}=\begin{bmatrix}\mathbf{0}&\tilde{\Pi}_{4}\\\ \mathbf{0}&-{\sf
diag}(\tilde{s}_{4,2})\end{bmatrix}$

Then it follows that
$\left(\tilde{\Pi}_{3},\tilde{s}_{3,1},\tilde{s}_{3,2},\tilde{t}_{3,1},\tilde{t}_{3,2}\right)$
is a feasible solution of Formulation 3. Finally we have:

$\displaystyle\langle
C_{aug},\tilde{\Pi}_{3}\rangle+\lambda\left[\|\tilde{s}_{3,1}\|_{1}+\|\tilde{s}_{3,2}\|_{1}+\|\tilde{t}_{3,1}\|_{1}+\|\tilde{t}_{3,2}\|_{1}\right]$
$\displaystyle=\langle
C_{aug},\tilde{\Pi}_{3}\rangle+\lambda\left[\|\tilde{s}_{4,1}\|_{1}+\|\tilde{s}_{4,2}\|_{1}\right]$
$\displaystyle=\langle
C,\tilde{\Pi}_{4}\rangle+\lambda\left[\|\tilde{s}_{4,1}\|_{1}+\|\tilde{s}_{4,2}\|_{1}\right]$
$\displaystyle<\langle
C,\Pi^{*}_{3,12}\rangle+\lambda\left[\|s^{*}_{3,1}-s^{*}_{3,2}\|_{1}+\|t^{*}_{3,1}-t^{*}_{3,2}\|_{1}\right]$
$\displaystyle=\langle
C_{aug},\Pi^{*}_{3}\rangle+\lambda\left[\|s^{*}_{3,1}\|_{1}+\|s^{*}_{3,2}\|_{1}+\|t^{*}_{3,1}\|_{1}+\|t^{*}_{3,2}\|_{1}\right]$

This contradicts the optimality of
$(\Pi^{*}_{3},s^{*}_{3,1},s^{*}_{3,2},t^{*}_{3,1},t^{*}_{3,2})$. This
completes the proof. ∎

###  A.3 Proof of continuous version

###### Proof.

In this proof we denote by $F_{1}$ the optimization problem of equation 2.3
and by $F_{2}$ the optimization problem equation 2.5. Let $\mu,\nu$ be two
absolutely continuous measures on $\mathbb{R}^{d}$. Moreover, we assume
$c(x,y)=\|x-y\|$ for some norm $\|\cdot\|$ on $\mathbb{R}^{d}$. We assume that
$\int\|x\|\nu(\mathrm{d}x),\int\|x\|\mu(\mathrm{d}x)<\infty$.

Step 1: Let $K_{\epsilon}$ be a compact set such that
$\int_{K_{\epsilon}}\|x\|\mu(\mathrm{d}x),\int_{K_{\epsilon}}\|x\|\nu(\mathrm{d}x)>1-\epsilon$.

Also, let $\tilde{K}_{\epsilon}=\\{x_{1},\dots,x_{n_{\epsilon}}\\}$ be a
maximal $\epsilon$-packing set of $K_{\epsilon}$. Starting from
$\tilde{K}_{\epsilon}$, define $\\{S_{1},\dots,S_{n_{\epsilon}}\\}$ as a
mutually disjoint covering of $K_{\epsilon}$ with internal points
$x_{1},\dots,x_{n_{\epsilon}}$ respectively, so that Diam$(S_{i})\leq
2\epsilon$. With $p_{i}=\int_{S_{i}}\mu(\mathrm{d}x)$,
$q_{i}=\int_{S_{i}}\nu(\mathrm{d}x)$ for $i=1,\dots,n_{\epsilon}$,
$p_{0}=\int_{K_{\epsilon}^{C}}\mu(\mathrm{d}x)$,
$q_{0}=\int_{K_{\epsilon}^{C}}\nu(\mathrm{d}x)$ and
$x_{0}=0\in\mathbb{R}^{d}$, define

$\displaystyle\mu_{\epsilon}$ $\displaystyle=$
$\displaystyle\sum_{0}^{n_{\epsilon}}p_{i}\delta_{x_{i}}$
$\displaystyle\nu_{\epsilon}$ $\displaystyle=$
$\displaystyle\sum_{0}^{n_{\epsilon}}q_{i}\delta_{x_{i}}$

A coupling $Q$ between two probability distributions is a joint distribution
with marginals as the given two distributions. The Wasserstein distance
between two distributions $P_{1}$ and $P_{2}$ is defined as:

$\displaystyle W_{1}(P_{1},P_{2})=\inf_{Q\in\mathscr{Q}(P_{1},P_{2})}\int
Q(x,y)\|x-y\|\mathrm{d}x\mathrm{d}y,$ (A.2)

where $\mathscr{Q}(P_{1},P_{2})$ is the collection of all couplings of $P_{1}$
and $P_{2}$.

Define $Q(x,y)=(\mathds{1}_{x=x_{0},y\in
K_{\epsilon}^{C}}+\sum_{i=1}^{n_{\epsilon}}\mathds{1}_{x=x_{i},y\in
S_{i}})\mu(\mathrm{d}y)$. Then $Q$ is a coupling between $\mu$ and
$\mu_{\epsilon}$. Therefore, clearly,

$\displaystyle
W_{1}(\mu,\mu_{\epsilon})\leq\int_{K_{\epsilon}^{C}}\|x\|\mu(\mathrm{d}x)+2\epsilon\left(\sum_{i=1}^{n_{\epsilon}}p_{i}\right)\leq
3\epsilon$ (A.3)

Similarly, $W_{1}(\nu,\nu_{\epsilon})\leq 3\epsilon$. Therefore
$\lim_{\epsilon\to 0}W_{1}(\nu,\nu_{\epsilon})=0$.

Moreover, $W_{1}(\mu,\nu)=\lim_{\epsilon\to
0}W_{1}(\mu_{\epsilon},\nu_{\epsilon})$, as
$W_{1}(\mu_{\epsilon},\nu_{\epsilon})-6\epsilon\leq W_{1}(\mu,\nu)\leq
W_{1}(\mu_{\epsilon},\nu_{\epsilon})+6\epsilon$ by triangle inequality.

Step 2: Let $S$ be an arbitrary measure with $\|S\|_{\mathrm{TV}}=2\gamma$, so
that $\mu+S$ is a probability measure with
$\int\|x\|(\mu+S)(\mathrm{d}x)<\infty$. Also, let us define
$\epsilon_{n}=2^{-(n+1)}$.

Let $S=S^{+}-S^{-}$, where $S^{+}$ and $S^{-}$ are positive measures on
$\mathbb{R}^{d}$. Then,
$\|S^{-}\|_{\mathrm{TV}}=\|S^{+}\|_{\mathrm{TV}}=\gamma$.

Clearly $(\mu-S^{-})/(1-\gamma),\mu,\nu,S^{+}/\gamma$ are tight probability
measures. So we can construct compact sets $K_{\epsilon_{n}}^{(1)}$, similar
to Step 1 to approximate all the four measures. Without loss of generality we
assume that $0\in K_{\epsilon_{n}}^{(1)}$ for all $n$. Moreover, we can also
construct approximate measures
$(\mu-S^{-})_{n}=((\mu-S^{-})/(1-\gamma))_{\epsilon_{n}}$ and
$(S^{+})_{n}=(S^{+}/\gamma)_{\epsilon_{n}}$ defined as in Step 1.
$\mu_{n}=\mu_{\epsilon_{n}},\nu_{n}=\nu_{\epsilon_{n}}$ are defined similarly.
All four of the measures have support points in $K_{\epsilon_{n}}^{(1)}$.

Next, we define $(\mu+S)_{n}=\gamma(S^{+})_{n}+(1-\gamma)(\mu-S^{-})_{n}$.
Then by the construction, from (Villani, 2009),
$\lim_{n\to\infty}W_{1}((\mu+S)_{n},\mu+S)\to 0$ and thus
$\lim_{n\to\infty}W_{1}((\mu+S)_{n},\nu_{n})\to W_{1}(\mu+S,\nu)$. Therefore
we can define a signed measure $S_{n}=(\mu+S)_{n}-\mu_{n}$. Moreover,

$\displaystyle\|S_{n}\|_{\mathrm{TV}}$
$\displaystyle\leq\gamma\|(S^{+})_{n}\|_{\mathrm{TV}}+\|(1-\gamma)(\mu-S^{-})_{n}-\mu_{n}\|_{\mathrm{TV}}$
(A.4) $\displaystyle=2\gamma=\|S\|_{\mathrm{TV}}$ (A.5)

Note that $\mu_{n},\nu_{n},(\mu+S)_{n}$ put masses (sometimes zero masses) on
a common set of support points given by $\tilde{K}_{\epsilon_{n}}^{(1)}\subset
K_{\epsilon_{n}}^{(1)}$.

The $\tilde{K}_{\epsilon_{n}}^{(1)}$ is sequentially defined so that
$\tilde{K}_{\epsilon_{n+1}}^{(1)}$ is a refinement of
$\tilde{K}_{\epsilon_{n}}^{(1)}$. This can easily be achieved by the choice of
$\epsilon_{n}$ defined.

Consider $\tilde{s}_{n},\Pi_{n}$ such that

$\displaystyle
F_{1}(\mu_{n},\nu_{n})=\int\|x-y\|\Pi_{n}(\mathrm{d}x\mathrm{d}y)+\lambda\|\tilde{s}_{n}\|_{\mathrm{TV}}$
(A.6)

By the discrete nature of $\mu_{n},\nu_{n}$, using the proof of the discrete
part $F_{1}(\mu_{n},\nu_{n})=F_{2}(\mu_{n},\nu_{n})$.

Since, $\min\\{\|x-y\|,2\lambda\\}$ is a metric, whenever $\|x-y\|$ is,
therefore, it is easy to check that
$F_{2}(\mu,\nu)=\lim_{n}F_{2}(\mu_{n},\nu_{n})=\lim_{n}F_{1}(\mu_{n},\nu_{n})$.

Moreover, by construction,
$F_{1}(\mu_{n},\nu_{n})\leq\int\|x-y\|\Pi(\mathrm{d}x\mathrm{d}y)+\lambda\|S_{n}\|_{\mathrm{TV}}$
for any arbitrary coupling $\Pi$ of $\mu$ and $\mu+S$. Also
$\lim_{n}W_{1}(\mu_{n},\mu),W_{1}(\mu+S,(\mu+S)_{n})\to 0$.

Thus, combining the above result with equation A.4, we get

$\displaystyle\lim_{n}F_{1}(\mu_{n},\nu_{n})\leq\int\|x-y\|\tilde{\Pi}(\mathrm{d}x\mathrm{d}y)+\lambda\|S\|_{\mathrm{TV}}$

for any coupling $\tilde{\Pi}$ of $\mu$ and $\mu+S$.

Therefore, $F_{2}(\mu,\nu)\leq F_{1}(\mu,\nu)$.

Step 3: Consider $\tilde{s}_{n}$ defined in equation A.6. As $\tilde{s}_{n}$
has support in the compact sets $K_{\epsilon_{n}}^{(1)}$ defined in Step 2,
therefore, $\\{\mu_{n}+\tilde{s}_{n}\\}_{n\geq 1}$ are tight measures.

Therefore, by Prokhorov’s Theorem for equivalence of sequential compactness
and tightness for a collection of measures, there exists a probability measure
$\mu\oplus s$ and a subsequence $\\{n_{k}\\}_{k\geq 1}$ such that
$\mu_{n_{k}}+\tilde{s}_{n_{k}}$ converges weakly to $\mu\oplus s$. Moreover,
by construction
$\lim_{R\to\infty}\limsup_{n\to\infty}\bigintss_{\|x\|>R}\|x\|(\mu_{n}+\nu_{n})(\text{d}x)=0$
and so
$\lim_{R\to\infty}\limsup_{n\to\infty}\bigintss_{\|x\|>R}\|x\|(\mu_{n}+\tilde{s}_{n})(\text{d}x)=0$.

Thus, by Definition 6.8 part (iii) and Theorem 6.9 of (Villani, 2009),
$W_{1}(\mu_{n_{k}}+\tilde{s}_{n_{k}},\mu\oplus s)\to 0$. Moreover,
$W_{1}(\mu_{n_{k}},\mu)\to 0$. Therefore
$\|\tilde{s}_{n_{k}}\|_{\mathrm{TV}}\to\|\mu\oplus s-\mu\|_{\mathrm{TV}}$.
Thus,
$W_{1}(\mu_{n_{k}}+\tilde{s}_{n_{k}},\nu_{n_{k}})+\lambda\|\tilde{s}_{n_{k}}\|_{\mathrm{TV}}\to
W_{1}(\mu\oplus s,\nu)+\lambda\|\mu\oplus s-\mu\|_{\mathrm{TV}}$. But by the
proof of the discrete part,
$W_{1}(\mu_{n_{k}}+\tilde{s}_{n_{k}},\nu_{n_{k}})+\lambda\|\tilde{s}_{n_{k}}\|_{\mathrm{TV}}=F_{1}(\mu_{n_{k}},\nu_{n_{k}})=F_{2}(\mu_{n_{k}},\nu_{n_{k}})\to
F_{2}(\mu,\nu)$. Therefore, with $s=\mu\oplus s-\mu$,
$W_{1}(\mu+s,\nu)+\lambda\|s\|_{\mathrm{TV}}=F_{2}(\mu,\nu)$.

Therefore, $F_{2}(\mu,\nu)=\limsup_{n\to\infty}F_{1}(\mu_{n},\nu_{n})\geq
F_{1}(\mu,\nu)$. Thus the equality holds.

∎

##  Appendix B Proof of Theorem 2.1

###### Proof.

The proof is immediate from the Formulation 1. Recall that the Formulation 1
can restructured as:

$\mathrm{ROBOT}(\tilde{\mu},\nu)=\inf_{P}\left\\{\mathrm{OT}(P,\nu)+\lambda\|P-\tilde{\mu}\|_{\mathrm{TV}}\right\\}\,.$

where the infimum is taking over all measure dominated by some common measure
$\sigma$ (with respect to which $\mu,\mu_{c},\nu$ are dominated). Hence,

$\mathrm{ROBOT}(\tilde{\mu},\nu)\leq\mathrm{OT}(P,\nu)+\lambda\|P-\tilde{\mu}\|_{\mathrm{TV}}$

for any particular choice of $P$. Taking $P=\mu$ we get that

$\displaystyle\mathrm{ROBOT}(\tilde{\mu},\nu)$
$\displaystyle\leq\mathrm{OT}(\mu,\nu)+\lambda\|\mu-\tilde{\mu}\|_{\mathrm{TV}}$
$\displaystyle=\mathrm{OT}(\mu,\nu)+\lambda{\epsilon}\|\mu-\mu_{c}\|_{\mathrm{TV}}$

Taking $P=\nu$ we get
$\mathrm{ROBOT}(\tilde{\mu},\nu)\leq\lambda\|\nu-\tilde{\mu}\|_{\mathrm{TV}}$
and finally taking $P=\tilde{\mu}$ we get
$\mathrm{ROBOT}(\tilde{\mu},\nu)\leq\mathrm{OT}(\tilde{\mu},\nu)$. This
completes the proof. ∎

##  Appendix C Proof of Lemma 3.2

As defined in the main text, let $\Pi^{*}_{2}$ be the optimal solution of
equation 2.7 and $\Pi^{*}_{2,\alpha}$ be the optimal solution of equation 3.1.
Then by Proposition 4.1 from Peyré & Cuturi (2018) we conclude:

$\Pi^{*}_{2,\alpha}\overset{\alpha\to 0}{\longrightarrow}\Pi^{*}_{2}\,.$ (C.1)

Now we have defined
$\left(\Pi^{*}_{1,\alpha},\mathbf{s}^{*}_{1,\alpha}\right)$ as the
_approximate_ solution of equation 2.6 obtained via Algorithm 1 from
$\Pi^{*}_{2,\alpha}$. Note that we can think of Algorithm 1 as a map from
${\mathbb{R}}^{m\times n}$ to
${\mathbb{R}}^{(m+n)\times(m+n)}\times{\mathbb{R}}^{m}$. Define this map as
$F$.

$F(\Pi_{2})\mapsto(\Pi_{1},\mathbf{s}_{1})$

Hence, by our notation,
$\left(\Pi^{*}_{1,\alpha},\mathbf{s}^{*}_{1,\alpha}\right)=F(\Pi^{*}_{2,\alpha})$
and $\left(\Pi^{*}_{1},\mathbf{s}^{*}_{1}\right)=F(\Pi^{*}_{2})$. Now if we
show that $F$ is a continuous map, then by continuous mapping theorem, it is
also immediate from equation C.1 that:

$\displaystyle F(\Pi^{*}_{2,\alpha})\overset{\alpha\to
0}{\longrightarrow}F(\Pi^{*}_{2})\,.$

which implies:

$\displaystyle\Pi^{*}_{1,\alpha}$ $\displaystyle\overset{\alpha\to
0}{\longrightarrow}\Pi^{*}_{1}$ $\displaystyle\mathbf{s}^{*}_{1,\alpha}$
$\displaystyle\overset{\alpha\to 0}{\longrightarrow}\mathbf{s}^{*}_{1}\,.$

which will complete the proof. Therefore all we need to show is that $F$ is a
continuous map. Towards that direction, first fix a sequence of matrices
$\\{\bar{\Pi}_{2,i}\\}_{i\in\mathbb{N}}\to\bar{\Pi}_{2}$. Define
$F(\bar{\Pi}_{2,i})=\left(\bar{\Pi}_{1,i},\bar{\mathbf{s}}_{1,i}\right)$ and
$F(\bar{\Pi}_{2})=\left(\bar{\Pi}_{1},\bar{\mathbf{s}}_{1}\right)$. By Step 3
- Step 5 of Algorithm 1, we obtain $\bar{\Pi}_{1,i}$ by first setting
$\bar{\Pi}_{1,i,12}=\bar{\Pi}_{2,i}$ and for each of the columns of
$\bar{\Pi}_{1,i,12}$, dumping the sum of its entries for which the cost is
$>2\lambda$ to the diagonals of $\bar{\Pi}_{1,i,22}$. Also, we have all the
entries of the first $n$ columns of $\bar{\Pi}_{1,i}$ to be $0$. In step 6 of
Algorithm 1, we obtain $\mathbf{s}_{1,i}$ by taking the negative of the sum of
the elements of each rows of $\bar{\Pi}_{1,i,12}$ for which the cost is
$>2\lambda$. Note that these operations (Step 3 - Step 6 of Algorithm 1) are
continuous. Therefore we conclude:

  1. 1.

$0=\bar{\Pi}_{1,i,11}\to\bar{\Pi}_{1,11}=0\,.$

  2. 2.

$0=\bar{\Pi}_{1,i,21}\to\bar{\Pi}_{1,21}=0\,.$

  3. 3.

$\bar{\Pi}_{1,i,12}=\bar{\Pi}_{2,i}\odot\mathds{1}_{\mathcal{I}^{c}}\to\bar{\Pi}_{2}\odot\mathds{1}_{\mathcal{I}^{c}}=\bar{\Pi}_{1,12}\,.$

  4. 4.

$\displaystyle\bar{\Pi}_{1,i,22}$ $\displaystyle={\sf
diag}\left(\mathbf{1}^{\top}\left(\bar{\Pi}_{2,i}\odot\mathds{1}_{\mathcal{I}}\right)\right)$
$\displaystyle\to{\sf
diag}\left(\mathbf{1}^{\top}\left(\bar{\Pi}_{2}\odot\mathds{1}_{\mathcal{I}}\right)\right)$
$\displaystyle=\bar{\Pi}_{1,22}\,.$

  5. 5.

$\displaystyle\mathbf{s}_{1,i}$
$\displaystyle=-\left(\bar{\Pi}_{i,n}\odot\mathds{1}_{\mathcal{I}}\right)\mathbf{1}$
$\displaystyle\to-\left(\bar{\Pi}_{2}\odot\mathds{1}_{\mathcal{I}}\right)\mathbf{1}=\mathbf{s}_{1}\,.$

where $A\odot B$ denotes the Hadamard product (element-wise multiplication)
between two matrices. Hence we have established:

$\displaystyle F(\bar{\Pi}_{2,i})$
$\displaystyle=\left(\bar{\Pi}_{1,i},\bar{\mathbf{s}}_{1,i}\right)$
$\displaystyle\overset{n\to\infty}{\longrightarrow}\left(\bar{\Pi}_{1},\bar{\mathbf{s}}_{1}\right)$
$\displaystyle=F(\bar{\Pi}_{2})\,.$

This completes the proof of continuity of $F$.

##  Appendix D Proof of auxiliary lemmas

###  D.1 Proof of Lemma A.1

###### Proof.

The fact that $\Pi^{*}_{1,11}=\Pi^{*}_{1,21}=\mathbf{0}$ follows from the fact
that $\Pi^{*}_{1}\succeq 0$ and $\Pi^{*}_{1}\mathbf{1}=\mathbf{Q}$. To prove
that $\Pi^{*}_{1,22}$ is diagonal, we use the fact that the any diagonal entry
the cost matrix is $0$. Now suppose $\Pi^{*}_{1,22}$ is not diagonal. Then
define a matrix $\hat{\Pi}$ as following: set
$\hat{\Pi}_{11}=\hat{\Pi}_{21}=\mathbf{0}$, $\hat{\Pi}_{12}=\Pi^{*}_{1,12}$
and:

$\hat{\Pi}_{22}(i,j)=\begin{cases}\sum_{k=1}^{m}\Pi^{*}_{1,22}(k,i),&\text{if
}j=i\\\ 0,&\text{if }j\neq i\end{cases}$

Also define $\hat{s}=s^{*}_{1}$ and $\hat{t}$ as
$\hat{t}(i)=\hat{\Pi}_{22}(i,i)$. Then clearly $(\hat{\Pi},\hat{s},\hat{t})$
is a feasible solution of Formulation 1. Note that:

$\|\hat{t}\|_{1}=1^{\top}\hat{\Pi}_{22}1=1^{\top}\Pi^{*}_{1,22}1=\|t^{*}_{1}\|_{1}$

and by our construction $\langle C_{aug},\hat{\Pi}\rangle<\langle
C_{aug},\Pi^{*}_{1}\rangle$. Hence $(\hat{\Pi},\hat{s},\hat{t})$ reduces the
value of the objective function of Formulation 1 which is a contradiction.
This completes the proof. ∎

###  D.2 Proof of Lemma A.2

###### Proof.

  1. 1.

Suppose $\Pi^{*}_{1}(i,j)>0$. Then dump this mass to $s^{*}_{1}(j)$ and make
it $0$. In this way $\langle C_{aug},\Pi^{*}_{1}\rangle$ will decrease by
$>2\lambda\Pi^{*}_{1}(i,j)$ and the regularizer value will increase by atmost
$2\lambda\Pi^{*}_{1}(i,j)$, resulting in overall reduction in the objective
value, which leads to a contradiction.

  2. 2.

Suppose each entry of $i^{th}$ row of $C$ is $<2\lambda$. Then if
$s^{*}_{1}(i)>0$, we can distribute this mass in the $i^{th}$ row such that,
$s^{*}_{1}(i)=a_{1}+a_{2}+\dots+a_{m}$ with the condition that
$t^{*}_{1}(j)\geq a_{j}$. Now we reduce $t^{*}_{1}$ as:

$t^{*}_{1}(j)\leftarrow t^{*}_{1}(j)-a_{j}$

Hence the value $\langle C_{aug},\Pi^{*}_{1}(i,j)\rangle$ will increase by a
value $<2\lambda s^{*}_{1}(i)$ but the value of regularizer will decrease by
the value of $2\lambda s^{*}_{1}(i)$, resulting in overall decrease in the
value of objective function.

  3. 3.

Same as proof of part (2) by interchanging row and column in the argument.

  4. 4.

Suppose not. Then choose ${\epsilon}<s^{*}_{1}(i)\wedge t^{*}_{1}(j)$, Add
${\epsilon}$ to $\Pi^{*}_{1}(i,j)$. Hence the cost function value $\langle
C_{aug},\Pi^{*}_{1}\rangle$ will increase by $<2\lambda{\epsilon}$ but the
regularizer value will decrease by $2\lambda{\epsilon}$, resulting in overall
decrease in the objective function.

∎

###  D.3 Proof of Lemma A.4

###### Proof.

For the notational simplicity, we drop the subscript $4$ now as we will only
deal with the solution of Formulation 4 and there will be no ambiguity. We
prove the Lemma by contradiction. Suppose $s^{*}_{1,i}>0$. Then we show one
can come up with another solution $(\tilde{\Pi},\tilde{s}_{1},\tilde{s}_{2})$
of Formulation 4 such that it has lower objective value. To construct this new
solution, make:

$\tilde{s}_{1,j}=\begin{cases}s^{*}_{1,j},&\text{if }j\neq i\\\ 0,&\text{if
}j=i\end{cases}$

Now to change the optimal transport plan, we will only change $i^{th}$ row of
$\Pi^{*}$. We subtract $a_{1},a_{2},\dots,a_{n}\geq 0$ from $i^{th}$ column of
$\Pi^{*}$ in such a way, such that none of the elements are negative. Hence
the column sum will be change, i.e. the value of $\tilde{s}_{2}$ will be:

$\tilde{s}_{2,j}=s^{*}_{2,j}-a_{j}\ \ \ \forall 1\leq j\leq n\,.$

Now clearly from our construction:

$\langle C,\tilde{\Pi}\rangle\leq\langle C,\Pi^{*}\rangle$

For the regularization part, note that, as we only reduced $i^{th}$ element of
$s^{*}_{1}$, we have $\|\tilde{s}_{1}\|_{1}=\|s^{*}_{1}\|_{1}-s^{*}_{1,i}$.
And by simple triangle inequality,

$\|\tilde{s}_{2}\|_{1}\leq\|s^{*}_{2}\|_{1}+\|a_{1}\|_{1}=\|s^{*}_{2}\|_{1}+s^{*}_{1,i}$

by construction $a_{i}$’s, as $a_{i}\geq 0$ and $\sum_{i}a_{i}=s^{*}_{1,i}$.
Hence we have:

$\|\tilde{s}_{1}\|_{1}+\|\tilde{s}_{2}\|_{1}\leq\|s^{*}_{1}\|_{1}-s^{*}_{1,i}+\|s^{*}_{2}\|_{1}+s^{*}_{1,i}=\|s^{*}_{1}\|_{1}+\|s^{*}_{2}\|_{1}\,.$

Hence the value corresponding to regularizer will also decrease. This
completes the proof. ∎

###  D.4 Proof of Lemma A.6

###### Proof.

We prove this lemma by contradiction. Suppose $\Pi^{*}_{3}$ does not have the
structure mentioned in the statement of Lemma. Construct another transport
plan for Formulation 3 $\tilde{\Pi}_{3}$ as follows: Keep
$\tilde{\Pi}_{3,12}=\Pi^{*}_{3,12}$ and set $\tilde{\Pi}_{3,12}=\mathbf{0}$.
Construct the other parts as:

$\displaystyle\tilde{\Pi}_{3,11}(i,j)=$
$\displaystyle\begin{cases}\sum_{k=1}^{m}\Pi^{*}_{3,11}(i,k)+\sum_{k=1}^{n}\Pi^{*}_{3,21}(k,i),&\text{if
}i=j\\\ 0,&\text{if }i\neq j\end{cases}$

and

$\tilde{\Pi}_{3,22}(i,j)=\begin{cases}\sum_{k=1}^{n}\Pi^{*}_{3,22}(k,i),&\text{if
}i=j\\\ 0,&\text{if }i\neq j\end{cases}$

It is immediate from the construction that:

$\langle C_{aug},\tilde{\Pi}_{3}\rangle\leq\langle C_{aug},\Pi^{*}_{3}\rangle$

As for the regularization term: Note the by our construction $\tilde{s}_{4}$
will be same as $s_{4}^{*}$ as column sum of $\tilde{\Pi}_{3,22}$ is same as
$\Pi^{*}_{3,22}$. For the other three:

$\tilde{s}_{3}(i)=\tilde{\Pi}_{3,11}(i,i)=\sum_{k=1}^{m}\Pi^{*}_{3,11}(i,k)+\sum_{k=1}^{n}\Pi^{*}_{3,21}(k,i)$
$\tilde{s}_{2}(i)=\tilde{\Pi}_{3,22}(i,i)=\sum_{k=1}^{n}\Pi^{*}_{3,22}(k,i)$

and hence by construction:

$\|\tilde{s}_{2}\|_{1}=\mathbf{1}^{\top}\Pi^{*}_{3,22}\mathbf{1}=\|s^{*}_{2}\|_{1}-\mathbf{1}^{\top}\Pi^{*}_{3,21}\mathbf{1}\,.$
$\|\tilde{s}_{3}\|_{1}=\mathbf{1}^{\top}\Pi^{*}_{3,11}\mathbf{1}+\mathbf{1}^{\top}\Pi^{*}_{3,21}\mathbf{1}=\|s^{*}_{3}\|_{1}$

And also by our construction, $\tilde{s}_{1}=s^{*}_{1}+c$ where
$c=(\Pi^{*}_{3,21})^{\top}\mathbf{1}$. As a consequence we have
$\|c\|_{1}=\mathbf{1}^{\top}\Pi^{*}_{3,21}\mathbf{1}$. Then it follows:

$\displaystyle\sum_{i=1}^{4}\|\tilde{s}_{i}\|_{1}$
$\displaystyle=\|s^{*}_{1}+c\|+\|s^{*}_{2}\|_{1}-\mathbf{1}^{\top}\Pi^{*}_{3,21}\mathbf{1}+\|s^{*}_{3}\|_{1}+\|s^{*}_{4}\|_{1}$
$\displaystyle\leq\sum_{i=1}^{4}\|s^{*}_{i}\|_{1}+\|c\|_{1}-\mathbf{1}^{\top}\Pi^{*}_{3,21}\mathbf{1}$
$\displaystyle=\sum_{i=1}^{4}\|s^{*}_{i}\|_{1}$

So the objective value is overall reduced. This contradicts the optimality of
$\Pi^{*}_{3}$ which completes the proof. ∎

##  Appendix E Change of support of outliers with respect to $\lambda$

For any $\lambda$, define the set
$\mathcal{I}_{\lambda}=\\{(i,j):C_{i,j}>2\lambda\\}$, i.e.
$\mathcal{I}_{\lambda}$ denotes the costs which exceeds the threshold
$2\lambda$. As before we define by $C_{\lambda}$ to be truncated cost $C\wedge
2\lambda$. Denote by $\pi_{\lambda}$ to be the optimal transport plan with
respect to $C_{\lambda}$ and the marginal measures $\mu,\nu$. Borrowing our
notations from previous theorems, we define a "slack vector"
$\mathbf{s}_{\lambda}$ as:

$\mathbf{s}_{\lambda}(i)=\sum_{j=1}^{n}\pi_{\lambda}(i,j)\mathds{1}_{C(i,j)>2\lambda}=\sum_{j:(i,j)\in\mathcal{I}_{\lambda}}\pi_{\lambda}(i,j)\,.$

And we define the observation $i_{0}$ to be an outlier if
$\mathbf{s}_{\lambda}(i_{0})>0$. It is immediate that for any
$\lambda_{1}<\lambda_{2}$,
$\mathcal{I}_{\lambda_{1}}\supseteq\mathcal{I}_{\lambda_{2}}$. We goal is to
establish the following theorem:

######  Theorem E.1.

For any $\lambda_{1}<\lambda_{2}$, if $\mathbf{s}_{\lambda_{2}}(i_{0})>0$,
then $\mathbf{s}_{\lambda_{1}}(i_{0})>0$, i.e. if a point is selected as
outlier for larger $\lambda$, then it is also selected as outlier for smaller
$\lambda$.

###### Proof.

Fix $\lambda_{1}<\lambda_{2}$. Note that for any $\pi\in\Pi(\mu,\nu)$ we have:

$\displaystyle\langle C_{\lambda_{2}}-C_{\lambda_{1}},\pi\rangle$
$\displaystyle=\sum_{(i,j)\in
I_{\lambda_{1}}\cap\mathcal{I}_{\lambda_{2}}^{c}}(C(i,j)-2\lambda_{1})\pi(i,j)$
$\displaystyle\qquad\qquad+2(\lambda_{2}-\lambda_{1})\sum_{(i,j)\in\mathcal{I}_{\lambda_{2}}}\pi(i,j)$
$\displaystyle:=T_{1}(\pi)+T_{2}(\pi)\,.$

Now as $\pi_{\lambda_{2}}$ is optimal with respect to $C_{\lambda_{2}}$ and
$\pi_{\lambda_{1}}$ is optimal with respect to $C_{\lambda_{1}}$ we have:

$\displaystyle\langle
C_{\lambda_{1}},\pi_{\lambda_{2}}\rangle+T_{1}(\pi_{\lambda_{2}})+T_{2}(\pi_{\lambda_{2}})$
$\displaystyle=\langle C_{\lambda_{2}},\pi_{\lambda_{2}}\rangle$
$\displaystyle\leq\langle C_{\lambda_{2}},\pi_{\lambda_{1}}\rangle$
$\displaystyle=\langle
C_{\lambda_{1}},\pi_{\lambda_{1}}\rangle+T_{1}(\pi_{\lambda_{1}})+T_{2}(\pi_{\lambda_{1}})$
$\displaystyle\leq\langle
C_{\lambda_{1}},\pi_{\lambda_{2}}\rangle+T_{1}(\pi_{\lambda_{1}})+T_{2}(\pi_{\lambda_{1}})$

Therefore we have:

$T_{1}(\pi_{\lambda_{2}})+T_{2}(\pi_{\lambda_{2}})\leq
T_{1}(\pi_{\lambda_{1}})+T_{2}(\pi_{\lambda_{1}})\,.$ (E.1)

But this is not enough. Note that we can further decompose $T_{1}$ (and
similarly $T_{2}$) as:

$\displaystyle T_{1,i}(\pi)$ $\displaystyle=\sum_{j:(i,j)\in
I_{\lambda_{1}}\cap\mathcal{I}_{\lambda_{2}}^{c}}\left(C_{i,j}-2\lambda_{1}\right)\pi(i,j)$
$\displaystyle T_{2,i}(\pi)$
$\displaystyle=2(\lambda_{2}-\lambda_{1})\sum_{j:(i,j)\in\mathcal{I}_{\lambda_{2}}}\pi(i,j)\,.$

Hence we have:

$T_{1}(\pi)=\sum_{i=1}^{n}T_{1,i}(\pi),\ \ \
T_{2}(\pi)=\sum_{i=1}^{n}T_{2,i}(\pi)\,.$

In equation E.1 we have established that
$T_{1}(\pi_{\lambda_{2}})+T_{2}(\pi_{\lambda_{2}})\leq
T_{1}(\pi_{\lambda_{1}})+T_{2}(\pi_{\lambda_{1}})\,.$ In addition if we can
show that:

$T_{1,i_{0}}(\pi_{\lambda_{1}})+T_{2,i_{0}}(\pi_{\lambda_{1}})=0\implies
T_{2,i_{0}}(\pi_{\lambda_{1}})=0\,.$ (E.2)

holds for all $1\leq i\leq n$ then we are done. This is because, suppose
$\mathbf{s}_{\lambda_{1}}(i_{0})=0$, Then

$T_{1,i_{0}}(\pi_{\lambda_{1}})+T_{2,i_{0}}(\pi_{\lambda_{1}})=0\,.$

This in turn by equation E.2 implies

$T_{2,i_{0}}(\pi_{\lambda_{2}})=0\,,$

i.e. $\mathbf{s}_{\lambda_{2}}(i_{0})=0$.

######  Lemma E.2.

Suppose $C$ is a $2\times 2$ cost matrix with all unequal cost:

$C=\begin{bmatrix}C_{11}&C_{12}\\\ C_{21}&C_{22}\end{bmatrix}$

If $C_{22}>2\lambda_{2}$ and $C_{21}<2\lambda_{1}$, then the following two
inequalities won’t occur simultaneously:

$\displaystyle C^{\lambda_{2}}_{11}+C^{\lambda_{2}}_{22}$ $\displaystyle\leq
C^{\lambda_{2}}_{12}+C^{\lambda_{2}}_{21}\,,$ $\displaystyle
C^{\lambda_{1}}_{12}+C^{\lambda_{1}}_{21}$ $\displaystyle\leq
C^{\lambda_{1}}_{11}+C^{\lambda_{1}}_{22}\,.$

Now suppose $i_{2}$ is an outlier with respect to $\lambda_{2}$ but not with
respect to $\lambda_{1}$. Then there exists $j_{1}$ and $j_{2}$ ($j_{1}\neq
j_{2}$) such that:

$C_{i_{2},j_{1}}<2\lambda_{1},C_{i_{2},j_{2}}>2\lambda_{2}$

such that $\pi^{\lambda_{2}}_{i_{2},j_{2}}>0$ and
$\pi^{\lambda_{1}}_{i_{2},j_{1}}>0$.

#### Case 1:

Now assume that we can find $i_{1}\neq i_{2}$ such that
$\pi^{\lambda_{2}}_{i_{1},j_{1}}>0$ and $\pi^{\lambda_{1}}_{i_{1},j_{2}}>0$.

Then $(i_{1},j_{2}),(i_{2},j_{1})\in{\sf supp}(\pi^{\lambda_{1}})$ and
$(i_{1},j_{1}),(i_{2},j_{2})\in{\sf supp}(\pi^{\lambda_{2}})$. Hence from
c-cyclical monotonicity properties of the support of the optimal transport
plan we have for $\pi^{\lambda_{1}}$:

$C^{\lambda_{1}}_{i_{1},j_{2}}+C^{\lambda_{1}}_{i_{2},j_{1}}\leq
C^{\lambda_{1}}_{i_{1},j_{1}}+C^{\lambda_{1}}_{i_{2},j_{2}}\,,$

and for $\pi^{\lambda_{2}}$:

$C^{\lambda_{2}}_{i_{1},j_{1}}+C^{\lambda_{2}}_{i_{2},j_{2}}\leq
C^{\lambda_{2}}_{i_{1},j_{2}}+C^{\lambda_{2}}_{i_{2},j_{1}}\,.$

which is a contradiction from Lemma E.2. This completes the proof.

#### Case 2: Need to be proved

Now we need to consider the other case, there does not exist any row
$i_{1}\neq i_{2}$ such that both $\pi^{\lambda_{2}}_{i_{1},j_{1}}>0$ and
$\pi^{\lambda_{1}}_{i_{1},j_{2}}>0$ occur simultaneously. This means that the
columns $j_{1},j_{2}$ are orthogonal, i.e.
$\langle\pi^{\lambda_{2}}_{:,j_{1}},\pi^{\lambda_{1}}_{:,j_{2}}\rangle=0$. ∎

###  E.1 Proof of Lemma E.2

As $C_{21}<2\lambda_{1}$ and $C_{22}>2\lambda_{2}$ we can modify the
inequalities in Lemma E.2 as:

$\displaystyle C^{\lambda_{2}}_{11}+2\lambda_{2}$ $\displaystyle\leq
C^{\lambda_{2}}_{12}+C_{21}\,,$ (E.3) $\displaystyle
C^{\lambda_{1}}_{12}+C_{21}$ $\displaystyle\leq
C^{\lambda_{1}}_{11}+2\lambda_{1}\,.$ (E.4)

Now as we have assume $C_{21}<2\lambda_{1}$, from equation E.3 we obtain:

$\displaystyle
2\lambda_{1}>C^{\lambda_{2}}_{11}-C^{\lambda_{2}}_{12}+2\lambda_{2}$
$\displaystyle\iff$ $\displaystyle
2(\lambda_{2}-\lambda_{1})<C^{\lambda_{2}}_{12}-C^{\lambda_{2}}_{11}\,.$ (E.5)

Hence $C^{\lambda_{2}}_{12}-C^{\lambda_{2}}_{11}>0$, which implies
$C_{11}<C_{12},C_{11}<2\lambda_{2}$ and also both $C_{11}$ and $C_{12}$ can
not lie within $(2\lambda_{1},2\lambda_{2})$. We divide the rest of the proofs
into four small cases:

#### Case 1:

Assume $2\lambda_{1}<C_{11}<2\lambda_{2},C_{12}>2\lambda_{2}$. In this case
from equation E.3 we have:

$\displaystyle C_{11}+2\lambda_{2}$ $\displaystyle\leq 2\lambda_{2}+C_{21}$

i.e. $C_{21}\geq C_{11}$ which is not possible as $C_{11}>2\lambda_{1}$ and
$C_{21}<2\lambda_{1}$.

#### Case 2:

Assume $C_{11}<2\lambda_{1},C_{12}>2\lambda_{2}$. Then from equation E.4 we
have $C_{21}\leq C_{11}$ and from equation E.3 we have: $C_{11}\leq C_{21}$
which cannot occur simultaneously.

#### Case 3:

Assume $C_{11}<2\lambda_{1}$ and $2\lambda_{1}<C_{12}<2\lambda_{2}$. Then from
equation E.3 and equation E.4 we have respectively:

$\displaystyle C_{11}+2\lambda_{2}$ $\displaystyle\leq C_{12}+C_{21}\,,$
$\displaystyle 2\lambda_{1}+C_{21}$ $\displaystyle\leq C_{11}+2\lambda_{1}$

From the second inequality we have $C_{21}\leq C_{11}$, which putting back in
the first inequality yields:

$C_{11}+2\lambda_{2}\leq C_{12}+C_{11}\implies C_{12}\geq 2\lambda_{2}$

which is a contradiction.

#### Case 4:

Assume $C_{11}<2\lambda_{1}$ and $C_{12}<2\lambda_{1}$. This form equation E.3
yields:

$\displaystyle C_{11}+2\lambda_{2}$ $\displaystyle\leq C_{12}+C_{21}$
$\displaystyle\implies C_{21}$ $\displaystyle\geq
C_{11}-C_{12}+2\lambda_{2}\,.$ (E.6)

Also from equation E.4 we have:

$\displaystyle C_{12}+C_{21}$ $\displaystyle\leq C_{11}+2\lambda_{1}$
$\displaystyle\implies C_{21}$ $\displaystyle\leq
C_{11}-C_{12}+2\lambda_{1}\,.$ (E.7)

From equation E.6 and equation E.7 we have:

$C_{11}-C_{12}+2\lambda_{2}\leq C_{11}-C_{12}+2\lambda_{1}$

i.e. $\lambda_{2}\leq\lambda_{1}$ which is a contradiction. This completes the
proof.

Table 5: Robust mean estimation with GANs using different distribution divergences. True mean is $\eta_{0}=\mathbf{0}_{5}$; sample size $n=1000$; contamination proportion ${\epsilon}=0.2$. We report results over 30 experiment restarts. Contamination | JS Loss | SH Loss | ROBOT
---|---|---|---
$\text{Cauchy}(0.1\cdot\mathbf{1_{5}},I_{5})$ | 0.2 $\pm$ 0.06 |  0.17 $\pm$ 0.04 |  0.17 $\pm$ 0.05
$\text{Cauchy}(0.5\cdot\mathbf{1_{5}},I_{5})$ | 0.3 $\pm$ 0.07 | 0.26 $\pm$ 0.05 |  0.25 $\pm$ 0.05
$\text{Cauchy}(1\cdot\mathbf{1_{5}},I_{5})$ | 0.45 $\pm$ 0.14 | 0.37 $\pm$ 0.06 |  0.36 $\pm$ 0.07
$\text{Cauchy}(2\cdot\mathbf{1_{5}},I_{5})$ | 0.39 $\pm$ 0.3 | 0.26 $\pm$ 0.06 |  0.2 $\pm$ 0.07

##  Appendix F Robust mean experiment with Cauchy distribution

In this section we present our results corresponding to the robust mean
estimation with the generative distribution $g_{\theta}(x)=x+\theta$ where
$x\sim\text{Cauchy}(0,1)$. As in Subsection 4.1, we assume that we have
observation $\\{x_{1},\dots,x_{n}\\}$ from a contaminated distribution
$(1-{\epsilon})\ \text{Cauchy}(\eta_{0},1)+{\epsilon}\
\text{Cauchy}(\eta_{1},1)$. For our experiments we take
$\eta_{0}=\mathbf{0}_{5}$ and vary
$\eta_{1}\in\left\\{0.1\cdot\mathbf{1_{5}},0.5\cdot\mathbf{1_{5}},1\cdot\mathbf{1_{5}},2\cdot\mathbf{1_{5}}\right\\}$
along wth ${\epsilon}=0.2$. We compare our method with Wu et al. (2020) and
results are presented in Table 5.

## References

  * Alvarez-Melis & Jaakkola (2018) David Alvarez-Melis and Tommi S Jaakkola.  Gromov-Wasserstein alignment of word embedding spaces.  _arXiv:1809.00013_ , 2018.
  * Arjovsky et al. (2017) Martin Arjovsky, Soumith Chintala, and Léon Bottou.  Wasserstein GAN.  _arXiv:1701.07875 [cs, stat]_ , January 2017.
  * Balaji et al. (2020) Yogesh Balaji, Rama Chellappa, and Soheil Feizi.  Robust optimal transport with applications in generative modeling and domain adaptation.  _Advances in Neural Information Processing Systems_ , 33, 2020.
  * Bassetti et al. (2006) Federico Bassetti, Antonella Bodini, and Eugenio Regazzini.  On minimum Kantorovich distance estimators.  _Statistics & Probability Letters_, 76(12):1298–1302, 2006.
  * Breunig et al. (2000) Markus M Breunig, Hans-Peter Kriegel, Raymond T Ng, and Jörg Sander.  Lof: identifying density-based local outliers.  In _Proceedings of the 2000 ACM SIGMOD international conference on Management of data_ , pp. 93–104, 2000.
  * Caffarelli & McCann (2010) Luis A Caffarelli and Robert J McCann.  Free boundaries in optimal transport and Monge-Ampere obstacle problems.  _Annals of Mathematics_ , pp. 673–730, 2010.
  * Chao et al. (2018) Gao Chao, Yao Yuan, and Zhu Weizhi.  Robust estimation via generative adversarial networks.  In _International Conference on Learning Representations_ , 2018.
  * Chizat et al. (2018) Lenaic Chizat, Gabriel Peyré, Bernhard Schmitzer, and François-Xavier Vialard.  Scaling algorithms for unbalanced optimal transport problems.  _Mathematics of Computation_ , 87(314):2563–2609, 2018.
  * Chizat. (2017) Lenaïc Chizat.  Unbalanced optimal transport: Models, numerical methods, applications.  _Numerical Analysis [math.NA]. Université Paris sciences et lettres_ , 2017.
  * Courty et al. (2014) Nicolas Courty, Rémi Flamary, and Devis Tuia.  Domain adaptation with regularized optimal transport.  In _Joint European Conference on Machine Learning and Knowledge Discovery in Databases_ , pp. 274–289. Springer, 2014.
  * Courty et al. (2017) Nicolas Courty, Rémi Flamary, Amaury Habrard, and Alain Rakotomamonjy.  Joint distribution optimal transportation for domain adaptation.  In _Advances in Neural Information Processing Systems_ , pp. 3730–3739, 2017.
  * Cuturi (2013) Marco Cuturi.  Sinkhorn Distances: Lightspeed Computation of Optimal Transport.  In C. J. C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K. Q. Weinberger (eds.), _Advances in Neural Information Processing Systems 26_ , pp. 2292–2300. Curran Associates, Inc., 2013.
  * Diamond & Boyd (2016) Steven Diamond and Stephen Boyd.  CVXPY: A Python-embedded modeling language for convex optimization.  _Journal of Machine Learning Research_ , 17(83):1–5, 2016.
  * Figalli (2010) Alessio Figalli.  The optimal partial transport problem.  _Archive for Rational Mechanics and Analysis_ , 195(2):533–560, 2010.
  * Flamary & Courty (2017) Rémi Flamary and Nicolas Courty.  POT Python optimal transport library, 2017.  URL https://pythonot.github.io/.
  * Genevay et al. (2016) Aude Genevay, Marco Cuturi, Gabriel Peyré, and Francis Bach.  Stochastic optimization for large-scale optimal transport.  In _Advances in Neural Information Processing Systems_ , pp. 3440–3448, 2016.
  * Ho et al. (2017) Nhat Ho, XuanLong Nguyen, Mikhail Yurochkin, Hung Hai Bui, Viet Huynh, and Dinh Phung.  Multilevel clustering via Wasserstein means.  In _International Conference on Machine Learning_ , pp. 1501–1509, 2017.
  * Huang et al. (2016) Gao Huang, Chuan Guo, Matt J Kusner, Yu Sun, Fei Sha, and Kilian Q Weinberger.  Supervised word mover’s distance.  In _Advances in Neural Information Processing Systems_ , pp. 4862–4870, 2016.
  * Huber & Ronchetti (2009) Peter J. Huber and Elvezio Ronchetti.  _Robust Statistics_.  Wiley Series in Probability and Statistics. Wiley, Hoboken, N.J, 2nd ed edition, 2009.  ISBN 978-0-470-12990-6.
  * Kantorovich (1942) Leonid Vitalievich Kantorovich.  On the translocation of masses.  In _Dokl. Akad. Nauk. USSR (NS)_ , volume 37, pp. 199–201, 1942\.
  * Kusner et al. (2015) Matt Kusner, Yu Sun, Nicholas Kolkin, and Kilian Weinberger.  From Word Embeddings To Document Distances.  In _International Conference on Machine Learning_ , pp. 957–966, June 2015.
  * Liero et al. (2018) Matthias Liero, Alexander Mielke, and Giuseppe Savaré.  Optimal entropy-transport problems and a new Hellinger–Kantorovich distance between positive measures.  _Inventiones Mathematicae_ , 211(3):969–1117, 2018.
  * Liu et al. (2008) Fei Tony Liu, Kai Ming Ting, and Zhi-Hua Zhou.  Isolation forest.  In _2008 eighth ieee international conference on data mining_ , pp. 413–422. IEEE, 2008.
  * Monge (1781) Gaspard Monge.  _Mémoire sur la théorie des déblais et des remblais_.  De l’Imprimerie Royale, 1781.
  * Nowozin et al. (2016) Sebastian Nowozin, Botond Cseke, and Ryota Tomioka.  f-GAN: Training generative neural samplers using variational divergence minimization.  In _Advances in Neural Information Processing Systems_ , pp. 271–279, 2016.
  * Pedregosa et al. (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830, 2011\.
  * Pele & Werman (2009) Ofir Pele and Michael Werman.  Fast and robust earth mover’s distances.  In _2009 IEEE 12th International Conference on Computer Vision_ , pp. 460–467. IEEE, 2009.
  * Peyré & Cuturi (2018) Gabriel Peyré and Marco Cuturi.  Computational Optimal Transport.  _arXiv:1803.00567 [stat]_ , March 2018.
  * Piccoli & Rossi (2014) Benedetto Piccoli and Francesco Rossi.  Generalized Wasserstein distance and its application to transport equations with source.  _Archive for Rational Mechanics and Analysis_ , 211(1):335–358, 2014.
  * Rousseeuw & Driessen (1999) Peter J Rousseeuw and Katrien Van Driessen.  A fast algorithm for the minimum covariance determinant estimator.  _Technometrics_ , 41(3):212–223, 1999.
  * Schölkopf et al. (1999) Bernhard Schölkopf, Robert C Williamson, Alexander J Smola, John Shawe-Taylor, John C Platt, et al.  Support vector method for novelty detection.  In _NIPS_ , volume 12, pp. 582–588. Citeseer, 1999.
  * Seguy et al. (2018) Vivien Seguy, Bharath Bhushan Damodaran, Rémi Flamary, Nicolas Courty, Antoine Rolet, and Mathieu Blondel.  Large-Scale Optimal Transport and Mapping Estimation.  _arXiv:1711.02283 [stat]_ , February 2018.
  * Shen & Sanghavi (2019) Yanyao Shen and Sujay Sanghavi.  Learning with bad training data via iterative trimmed loss minimization.  In _International Conference on Machine Learning_ , pp. 5739–5748. PMLR, 2019.
  * Solomon et al. (2015) Justin Solomon, Fernando De Goes, Gabriel Peyré, Marco Cuturi, Adrian Butscher, Andy Nguyen, Tao Du, and Leonidas Guibas.  Convolutional Wasserstein distances: Efficient optimal transportation on geometric domains.  _ACM Transactions on Graphics (TOG)_ , 34(4):1–11, 2015.
  * Srivastava et al. (2018) Sanvesh Srivastava, Cheng Li, and David B. Dunson.  Scalable Bayes via Barycenter in Wasserstein Space.  _arXiv:1508.05880 [stat]_ , June 2018.
  * Staerman et al. (2020) Guillaume Staerman, Pierre Laforgue, Pavlo Mozharovskyi, and Florence d’Alché Buc.  When OT meets MOM: Robust estimation of Wasserstein distance.  _arXiv:2006.10325_ , 2020.
  * Tagasovska & Lopez-Paz (2019) Natasa Tagasovska and David Lopez-Paz.  Single-model uncertainties for deep learning.  _Advances in Neural Information Processing Systems_ , 32:6417–6428, 2019.
  * Tong et al. (2020) Alexander Tong, Guy Wolf, and Smita Krishnaswamyt.  Fixing bias in reconstruction-based anomaly detection with lipschitz discriminators.  In _2020 IEEE 30th International Workshop on Machine Learning for Signal Processing (MLSP)_ , pp. 1–6. IEEE, 2020.
  * Villani (2009) C. Villani.  _Optimal Transport: Old and New. Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathemtical Sciences]_.  Springer, Berlin, 2009.
  * Wu et al. (2020) Kaiwen Wu, Gavin Weiguang Ding, Ruitong Huang, and Yaoliang Yu.  On Minimax Optimality of GANs for Robust Mean Estimation.  In _International Conference on Artificial Intelligence and Statistics_ , pp. 4541–4551, June 2020.
  * Yurochkin et al. (2019) Mikhail Yurochkin, Sebastian Claici, Edward Chien, Farzaneh Mirzazadeh, and Justin Solomon.  Hierarchical Optimal Transport for Document Representation.  _arXiv:1906.10827 [cs, stat]_ , June 2019.