# Sliced Optimal Partial Transport

Yikun Bai<EMAIL_ADDRESS>Bernard Schmitzer
<EMAIL_ADDRESS>Mathew Thorpe
<EMAIL_ADDRESS>The Alan Turing Institute  Soheil Kolouri
<EMAIL_ADDRESS>
###### Abstract

Optimal transport (OT) has become exceedingly popular in machine learning,
data science, and computer vision. The core assumption in the OT problem is
the equal total amount of mass in source and target measures, which limits its
application. Optimal Partial Transport (OPT) is a recently proposed solution
to this limitation. Similar to the OT problem, the computation of OPT relies
on solving a linear programming problem (often in high dimensions), which can
become computationally prohibitive. In this paper, we propose an efficient
algorithm for calculating the OPT problem between two non-negative measures in
one dimension. Next, following the idea of sliced OT distances, we utilize
slicing to define the sliced OPT distance. Finally, we demonstrate the
computational and accuracy benefits of the sliced OPT-based method in various
numerical experiments. In particular, we show an application of our proposed
Sliced-OPT in noisy point cloud registration.

**footnotetext: These authors contributed equally to this work.

##  1 Introduction

The Optimal Transport (OT) problem studies how to find the most cost efficient
way to transport one probability measure to another, and it gives rise to
popular probability metrics like the Wasserstein distance. OT has attracted
abundant attention in data science, statistics, machine learning, signal
processing and computer vision [50, 23, 40, 32, 1, 25, 38, 52, 14, 13]

A core assumption in the OT problem is the equal total amount of mass in the
source and target measures (e.g., probability measures). Many practical
problems, however, deal with comparing non-negative measures with varying
total amount of mass, e.g., shape analysis [49, 10], domain adaptation [18],
color transfer [11]. In addition, OT distances are often not robust to
outliers and noise, as transporting outliers could be prohibitively expensive
and might compromise the distance estimation. To address these issues, many
variants of the OT problem have been recently proposed, for example the
optimal partial transport (OPT) problem [7, 20, 21], the Hellinger–Kantorovich
distance [10, 37], and Kantorovich–Rubinstein norm [26, 36]. These variants
were subsequently unified under the name “unbalanced optimal tansport” [12,
37].

The computational complexity of linear programming for balanced and partial OT
problems is often a bottleneck for solving large scale problems. Different
approaches have been developed to address this issue. For instance, by
entropic regularization the problem becomes strictly convex and can be solved
with the celebrated Sinkhorn–Knopp algorithm [15, 48] which has been extended
to the unbalanced setting [11]. This approach can still be computationally
expensive for small regularization levels. Other strategies exploit specific
properties of ground costs. For example, if the ground cost is determined by
the unique path on a tree, the problem can be efficiently solved in the
balanced [42, 34] and the unbalanced setting [47]. In particular, balanced
1-dimensional transport problems with convex ground costs can be solved by the
north-west corner rule, which essentially amounts to sorting the support
points the two input measures.

Base on this, another popular method is the sliced OT approach [45, 31, 5],
which assumes the ground cost is consistent with the Euclidean distance (in
1-dimensional space). Furthermore, it has been shown [31, 33, 47] that the OT
distance in Euclidean space can be approximated by the OT distance in
1-dimensional Euclidean space. Inspired by these works, in this paper, we
propose the sliced version of OPT and an efficient computational algorithm for
empirical distributions with uniform weights, i.e., measures of the form
$\sum_{i=1}^{n}\delta_{x_{i}}$, where $\delta_{x}$ is the Dirac measure. Our
contributions in this paper can be summarized as follows:

  * •

We propose a primal-dual algorithm for 1-dimensional OPT with a quadratic
worst case time complexity and linear or quadratic complexity in practice.

  * •

In $d$-dimensional space, we propose a the Sliced-OPT (SOPT) distance. Similar
to the sliced OT distance, we prove that it satisfies the metric axioms and
propose an computational method based on our 1-dimensional OPT problem solver.

  * •

We demonstrate an application of SOPT in point cloud registration by proposing
a SOPT variant of the iterative closest point (ICP) algorithm. Our approach is
robust against noise. Also, we apply SOPT to a color adaptation problem.

##  2 Related Work

Linear programming. In the discrete case the Kantorovich formulation [28] of
OT problem is a (high-dimensional) linear program [29]. As shown in [7], OPT
can be formulated as a balanced OT problem by introducing reservoir points,
thus it could also be solved by linear programming. However, time complexity
is prohibitive for large datasets.

Entropy Regularization. Entropic regularization approaches add the transport
plan’s entropy to the OT objective function and then apply the Sinkhorn-Knopp
algorithm [48, 15]. The algorithm can be extended to the large-scale
stochastic [24] and unbalanced setting [2, 11]. For moderate regularization
these algorithms converge fast, however there is a trade-off between accuracy
versus stability and convergence speed for small regularization.

Sliced OT. Sliced OT techniques [45, 31, 5, 39, 33] rely on the closed form
solution for the balanced OT map in 1-dimensional Euclidean settings, i.e.,
the increasing re-arrangement function given by the north-west corner rule.
The main idea behind these methods is to calculate the expected OT distance
between the 1-dimensional marginal distributions (i.e., slices) of two
$d$-dimensional distributions. The expectation is numerically approximated via
a Monte Carlo integration scheme. Other extensions of these distances include
the generalized and the max-sliced Wasserstein distances [30, 16]. In the
unbalanced setting, and for a particular case of OPT, [4] propose a fast
(primal) algorithm, which has quadratic worst case time complexity, and often
linear complexity in practice. In particular, Bonneel et al. [4] assume that
all the mass in the source measure must be transported to the target measure,
i.e., no mass destruction happens in the source measure.

Other computational methods. When the transportation cost is a metric, network
flow methods [26, 41] can be applied. For metrics on trees an efficient
algorithm based on dynamic programming with time complexity
$\mathcal{O}(n\log^{2}n)$ is proposed in [47]. However, in high dimensions
existence (and identification) of an appropriate metric tree remains
challenging.

##  3 Background of Optimal (Partial) Transport

We first review the preliminary concepts of the OT and the OPT problems. In
what follows, given $\Omega\subset\mathbb{R}^{r},p\geq 1$, we denote by
$\mathcal{P}(\Omega)$ the set of Borel probability measures and by
$\mathcal{P}_{p}(\Omega)$ the set of probability measures with finite $p$’th
moment defined on a metric space $(\Omega,d)$.

Optimal transport. Given $\mu,\nu\in\mathcal{P}(\Omega)$, and a lower semi-
continuous function $c:\Omega^{2}\to\mathbb{R}_{+}$, the OT problem between
$\mu$ and $\nu$ in the Kantorovich formulation [28], is defined as:

$\displaystyle\text{OT}(\mu,\nu):=\inf_{\gamma\in\Gamma(\mu,\nu)}\int_{\Omega^{2}}c(x,y)\,\mathrm{d}\gamma(x,y),$
(1)

where $\Gamma(\mu,\nu)$ is the set of all joint probability measures whose
marginal are $\mu$ and $\nu$. Mathematically, we denote as
$\pi_{1\\#}\gamma=\mu,\pi_{2\\#}\gamma=\nu$, where $\pi_{1},\pi_{2}$ are
canonical projection maps, and for any (measurable) function
$f:\Omega^{2}\to\Omega$, $f_{\\#}\gamma$ is the push-forward measure defined
as $f_{\\#}\gamma(A)=\gamma(f^{-1}(A))$ for any Borel set $A\subset\Omega$.
When $c(x,y)$ is a $p$-th power of a metric, the $p$-th root of the induced
optimal value is the Wasserstein distance, a metric in
$\mathcal{P}_{p}(\Omega)$.

Optimal Partial Transport. The OPT problem, in addition to mass
transportation, allows mass destruction on the source and mass creation on the
target. Here the mass destruction and creation penalty will be linear. Let
$\mathcal{M}_{+}(\Omega)$ denote the set of all positive Radon measures
defined on $\Omega$, suppose $\mu,\nu\in\mathcal{M}_{+}(\Omega)$, and
$\lambda_{1},\lambda_{2}\geq 0$, the OPT problem is:

$\displaystyle\mathrm{OPT}_{\lambda_{1},\lambda_{2}}(\mu,\nu):=\inf_{\begin{subarray}{c}\gamma\in\mathcal{M}_{+}(\Omega^{2})\\\
\pi_{1\\#}\gamma\leq\mu,\pi_{2\\#}\gamma\leq\nu\end{subarray}}\int
c(x,y)\,\mathrm{d}\gamma+\lambda_{1}(\mu(\Omega)-\pi_{1\\#}\gamma(\Omega))+\lambda_{2}(\nu(\Omega)-\pi_{2\\#}\gamma(\Omega))$
(2)

where the notation $\pi_{1\\#}\gamma\leq\mu$ denotes that for any Borel set
$A\subseteq\Omega$, $\pi_{1\\#}\gamma(A)\leq\mu(A)$, and we say
$\pi_{1\\#}\gamma$ is dominated by $\mu$, analogously for
$\pi_{2\\#}\gamma\leq\nu$; the notation $\mu(\Omega)$ denotes the total mass
of measure $\mu$. We denote the set of such $\gamma$ by
$\Gamma_{\leq}(\mu,\nu)$. When the transportation cost $c(x,y)$ is a metric,
and $\lambda_{1}=\lambda_{2}$, $\mathrm{OPT}(\cdot,\cdot)$ defines a metric on
$\mathcal{M}_{+}(\Omega)$ (see [12, Proposition 2.10], [43, Proposition 5],
[35, Section 2.1] and [8, Theorem 4]). For finite $\lambda_{1}$ and
$\lambda_{2}$ let $\lambda=\frac{\lambda_{1}+\lambda_{2}}{2}$ and define:

$\displaystyle\mathrm{OPT}_{\lambda}(\mu,\nu):=$
$\displaystyle\mathrm{OPT}_{\lambda,\lambda}(\mu,\nu)$ (3) $\displaystyle=$
$\displaystyle\mathrm{OPT}_{\lambda_{1},\lambda_{2}}(\mu,\nu)-K_{\lambda_{1},\lambda_{2}}(\mu,\nu)$

where,

$\displaystyle
K_{\lambda_{1},\lambda_{2}}(\mu,\nu)=\frac{\lambda_{1}-\lambda_{2}}{2}\mu(\Omega)+\frac{\lambda_{2}-\lambda_{1}}{2}\nu(\Omega).$

Since for fixed $\mu$ and $\nu$, $K_{\lambda_{1},\lambda_{2}}$ is a constant,
and without the loss of generality, in the rest of the paper we only consider
$\mathrm{OPT}_{\lambda}(\mu,\nu)$. The case where $\lambda_{i}$s are not
finite is discussed in Section 4.1.

Various equivalent formulations of OPT (3) have appeared in prior work, e.g.,
[21, 20, 43], which were later unified as a special case of unbalanced OT [12,
37]. We provide a short summary of these formulations and their relationship
with unbalanced OT in the appendix. OPT has several desirable theoretical
properties. For instance, by [43, Proposition 5], minimizing $\gamma$ exist
and they are concentrated on $c$-cyclical monotone sets. More concretely, we
have the following proposition.

######  Proposition 3.1.

Let $\gamma^{*}$ be a minimizer in (3), then the support of $\gamma^{*}$
satisfies the $c$-cyclical monotonicity property: for any $n\in\mathbb{N}$,
any $\\{(x_{i},y_{i})\\}_{i=1}^{n}\subset\mathrm{supp}(\gamma^{*})$ and any
permutation $\sigma:[1:n]\to[1:n]$ we have

$\sum_{i=1}^{n}c(x_{i},y_{i})\leq\sum_{i=1}^{n}c(x_{i},y_{\sigma(i)}).$

In particular, in one dimension, for $c(x,y)=f(x-y)$ where
$f:\mathbb{R}\to\mathbb{R}$ is a convex function, c-cyclical monotonicity is
equivalent to

$\displaystyle(x_{1},y_{1}),(x_{2},y_{2})\in\mathrm{supp}(\gamma^{*})\qquad\Rightarrow\qquad[x_{1}\leq
x_{2}\text{ and }y_{1}\leq y_{2}]\text{ or }[x_{1}\geq x_{2}\text{ and
}y_{1}\geq y_{2}].$

###### Proof.

Let $\hat{\gamma}$ be optimal for the extended balanced problem of appendix
section B, (11), and let $\gamma$ be the restriction of this measure to
$\Omega\times\Omega$. Since restriction preserves optimality [54, Theorem
4.6], $\gamma$ must be an optimal plan between $\pi_{1\\#}\gamma$ and
$\pi_{2\\#}\gamma$ with respect to the (non-extended) cost $c$ on
$\Omega\times\Omega$. Therefore, it must be supported on a $c$-cyclically
monotone set [54, Theorem 5.10]. In one dimension, for costs of the form
$c(x,y)=f(x-y)$ for convex $f$, $c$-cyclical monotonicity reduces to standard
monotonicity, see for instance [46, Theorem 2.9]. ∎

To further simplify (3), we show in the following lemma that the support of
the optimal $\gamma$ does not contain pairs of $(x,y)$ whose cost exceeds
$2\lambda$.

######  Lemma 3.2.

Let $\gamma^{*}$ be a minimizer in (3). Then, for all
$(x,y)\in\mathrm{supp}(\gamma^{*})$ we have $c(x,y)<2\lambda$.

###### Proof.

Pick $\gamma$ and let $S=\\{(x,y)\in\mathrm{supp}(\gamma):c(x,y)\geq
2\lambda\\}$. Define $\gamma^{\prime}$ as follows: for any Borel
$A\subset\Omega^{2}$, $\gamma^{\prime}(A)=\gamma(A\setminus S)$. Let

$\displaystyle C(\gamma)$ $\displaystyle:=\int
c(x,y)\,\mathrm{d}\gamma+\lambda\left[(\mu(\Omega)-\pi_{1\\#}\gamma(\Omega))+(\nu(\Omega)-\pi_{2\\#}\gamma(\Omega))\right]$
$\displaystyle=\int\left(c(x,y)-2\lambda\right)\,\mathrm{d}\gamma+\lambda(\mu(\Omega)+\nu(\Omega)),$

which is the objective function for the OPT problem defined in (3), and the
second line follows from the fact
$\gamma(\Omega^{2})=(\pi_{1})_{\\#}\gamma(\Omega)=(\pi_{2})_{\\#}\gamma(\Omega)$.
Then we have

$\displaystyle
C(\gamma)-C(\gamma^{\prime})=\int_{S}(c(x,y)-2\lambda)\,\mathrm{d}\gamma(x,y)\geq
0$

That is, for any $\gamma$ who has support on point $(x,y)$ such that
$c(x,y)\geq 2\lambda$, we can find a better transportation plan
$\gamma^{\prime}$. ∎

##  4 Empirical Optimal Partial Transport

In $\mathbb{R}^{r}$, suppose $\mu,\nu$ are $n$ and $m$-size empirical
distributions, i.e., $\mu=\sum_{i=1}^{n}\delta_{x_{i}}$ and
$\nu=\sum_{j=1}^{m}\delta_{y_{j}}$. The OPT problem (3), denoted as
$\mathrm{OPT}(\\{x_{i}\\}_{i=1}^{n},\\{y_{j}\\}_{j=1}^{m})$ can be written as

$\displaystyle\mathrm{OPT}(\\{x_{i}\\}_{i=1}^{n},\\{y_{j}\\}_{j=1}^{m}):=\min_{\gamma\in\Gamma_{\leq}(\mu,\nu)}\sum_{i,j}c(x_{i},y_{j})\gamma_{ij}+\lambda(n+m-2\sum_{i,j}\gamma_{ij})$
(4)

where

$\Gamma_{\leq}(\mu,\nu):=\\{\gamma\in\mathbb{R}_{+}^{n\times m}:\gamma
1_{m}\leq 1_{n},\gamma^{T}1_{n}\leq 1_{m}\\},$

and $1_{n}$ denotes the $n\times 1$ vector whose entries are $1$ and
analogously for $1_{m}$. We show that the optimal plan $\gamma$ for the
empirical OPT problem is induced by a 1-1 mapping. A similar result is known
for continuous measures $\mu$ and $\nu$, see [20, Proposition 2.4 and Theorem
2.6].

######  Theorem 4.1.

There exists an optimal plan for
$\mathrm{OPT}(\\{x_{i}\\}_{i=1}^{n},\\{y_{j}\\}_{j=1}^{m})$, which is induced
by a 1-1 mapping, i.e., $\gamma_{ij}\in\\{0,1\\},\forall i,j$ and in each row
and column, at most one entry of $\gamma$ is 1.

###### Proof.

We start by adapting the extension (11) to the concrete discrete setting
between empirical measures. Let $\hat{\Omega}=[1:m+n]$, and let
$\hat{c}:\hat{\Omega}\times\hat{\Omega}$ be given by

$\hat{c}(i,j)=\begin{cases}c(x_{i},y_{j})-2\lambda&\text{if }i\leq n,j\leq
m,\\\ 0&\text{otherwise.}\end{cases}$

Let $\hat{\mu}=\hat{\nu}=\sum_{i=1}^{m+n}\delta_{i}$. Then solving the
(balanced) optimal transport problem on $\hat{\Omega}$ between $\hat{\mu}$ and
$\hat{\nu}$ with respect to cost $\hat{c}$ is (as above) clearly equivalent to
the OPT problem, and an optimal OPT plan can be obtained by restricting an
optimal $\hat{\gamma}$ to the set $[1:n]\times[1:m]$. Note that here we have
simply split the mass on the isolated point $\hat{\infty}$ onto multiple
points, where each only carries unit mass. At the same time, the set
$\Gamma(\hat{\mu},\hat{\nu})$ are the doubly stochastic matrices and by the
Birkhoff-von-Neumann theorem its extremal points are permutation matrices.
Thus there always exists an optimal $\hat{\gamma}$ that is a permutation
matrix, and thus its restriction to $[1:n]\times[1:m]$ will only contain
entries $0$ or $1$, with at most one $1$ per row and column. ∎

Combining this theorem and the cyclic monotonicity of 1D OPT, we can restrict
the optimal mapping to strictly increasing maps.

######  Corollary 4.2.

For $\\{x_{i}\\}_{i=1}^{n},\\{y_{j}\\}_{j=1}^{m}$ sorted point lists in
$\mathbb{R}$, and a cost function $c(x,y)=h(x-y)$ where
$h:\mathbb{R}\to\mathbb{R}$ is strictly convex, the empirical OPT problem
$\mathrm{OPT}(\\{x_{i}\\}_{i=1}^{n},\\{y_{j}\\}_{j=1}^{m})$ can be further
simplified to

$\displaystyle\mathrm{OPT}(\\{x_{i}\\}_{i=1}^{n},\\{y_{j}\\}_{j=1}^{m}):=\min_{L}C(L)$
(5)

where

$\displaystyle
C(L):=\sum_{i\in\mathrm{dom}(L)}c(x_{i},y_{L[i]})+\lambda\left(n+m-2\left|\mathrm{dom}(L)\right|\right),$

and $L:[1:n]\to\\{-1\\}\cup[1:m]$, $\mathrm{dom}(L):=\\{i:L[i]\neq-1\\}$,
$L_{\mid\mathrm{dom}(L)}:[1:n]\hookrightarrow[1:m]$ is a strictly increasing
mapping.***Here, mapping a point to $\\{-1\\}$ corresponds to destroying the
point in the source, while unmatched points in the target are created. Hence,
$L$ uniquely represents the partial transport.

For convenience, we call any mapping $L:[1:n]\to\\{-1\\}\cup[1:m]$ a
“transportation plan” (or “plan” for short). Furthermore, since $L$ can be
represented by a vector, we do not distinguish $L(i)$, $L[i]$ and $L_{i}$. Of
course, there is a bijection between admissible $L$ and $\gamma$ and we use
this equivalence implicitly in the following.

Importantly, the OPT problem is a convex optimisation problem and therefore
has a dual form which is given by the following proposition.

######  Proposition 4.3.

The primal problem (4) admits the dual form

$\sup_{\begin{subarray}{c}\Phi\in\mathbb{R}^{n},\Psi\in\mathbb{R}^{m}\\\
\Phi_{i}+\Psi_{j}\leq c(x_{i},y_{j})\,\forall
i,j\end{subarray}}\sum_{i=1}^{n}\min\\{\Phi_{i},\lambda\\}+\sum_{j=1}^{m}\min\\{\Psi_{j},\lambda\\}.$

Moreover, the following are necessary and sufficient conditions for
$\gamma\in\Gamma_{\leq}(\mu,\nu)$, $\Phi\in\mathbb{R}^{n}$ and
$\Psi\in\mathbb{R}^{m}$ to be optimal for the primal and dual problems:

$\Phi_{i}+\Psi_{j}=c(x_{i},y_{j}),\forall(x_{i},y_{j})\in\mathrm{supp}(\gamma)$
---
$\Phi_{i}<\lambda\Rightarrow[\pi_{1\\#}\gamma]_{i}=1$ |  $\Psi_{j}<\lambda\Rightarrow[\pi_{2\\#}\gamma]_{j}=1$
$\Phi_{i}=\lambda\Rightarrow[\pi_{1\\#}\gamma]_{i}\in[0,1]$ |  $\Psi_{j}=\lambda\Rightarrow[\pi_{2\\#}\gamma]_{j}\in[0,1]$
$\Phi_{i}>\lambda\Rightarrow[\pi_{1\\#}\gamma]_{i}=0$ |  $\Psi_{i}>\lambda\Rightarrow[\pi_{2\\#}\gamma]_{j}=0$.

###### Proof.

For $\lambda_{1}=\lambda_{2}$ we can write (2) as

$\mathrm{OPT}_{\lambda}(\mu,\nu)=\mathrm{ET}(\mu,\nu;\mathcal{F},\mathcal{F})=\inf_{\gamma\geq
0}\int_{\Omega^{2}}c\mathrm{d}\gamma+F(\pi_{1\\#}\gamma|\mu)+F(\pi_{2\\#}\gamma|\nu)$

where $\mathcal{F}(\hat{\mu}\parallel\mu)$ is called Csiszàr $f-$divergence,
defined as

$\mathcal{F}(\hat{\mu}\parallel\mu):=\int
F(\sigma)d\mu+F^{\prime}_{\infty}\gamma^{\perp}(\Omega)=\begin{cases}\lambda\cdot(\mu(\Omega)-\hat{\mu}(\Omega))&\text{if
}0\leq\hat{\mu}\leq\mu\\\ +\infty&\text{otherwise}\end{cases},$

with the integrand $F$

$F(s)=\begin{cases}\lambda(1-s)&\text{if }s\in[0,1]\\\
+\infty&\text{else}\end{cases},$

and $\sigma,\mu^{\perp}$ is defined by Lebesgue’s decomposition theorem
$\mu=\frac{d\hat{\mu}}{d\mu}+\mu^{\perp}$;
$F^{\prime}_{\infty}:=\lim_{s\to\infty}\frac{F(s)}{s}$ is called recession
constant of $F$ (We refer [37, section 2.1] for more details).

By [37, Theorem 4.11] the dual of $\mathrm{ET}$ is

$\sup_{\begin{subarray}{c}\Phi\in\mathrm{L}^{1}(\mu),\Psi\in\mathrm{L}^{1}(\nu)\\\
\Phi\oplus\Psi\leq c\\\ \Phi,\Psi\text{ lsc and
bounded}\end{subarray}}-\int_{\Omega}F^{*}(-\Phi)\,\mathrm{d}\mu-\int_{\Omega}F^{*}(-\Psi)\,\mathrm{d}\nu$

where $F^{*}:\mathbb{R}\to(-\infty,\infty]$, called Legendre conjugate
function, is defined as

$F^{*}(r)=\sup_{s}(rs-F(s))=\max\\{-\lambda,r\\}.$

By [37, Theorem 4.6] the optimality conditions are

$\displaystyle\Phi\oplus\Psi$ $\displaystyle=c$
$\displaystyle\gamma\text{-a.e.}$ $\displaystyle-\Phi$
$\displaystyle\in\partial
F\left(\frac{\mathrm{d}\gamma_{1}}{\mathrm{d}\mu}\right),\quad\gamma_{1}=\pi_{1\\#}\gamma$
$\displaystyle\mu\text{-a.e.}$ (6) $\displaystyle-\Psi$
$\displaystyle\in\partial
F\left(\frac{\mathrm{d}\gamma_{2}}{\mathrm{d}\nu}\right),\quad\gamma_{2}=\pi_{2\\#}\gamma$
$\displaystyle\nu\text{-a.e.}$

We have

$\partial F(s)=\left\\{\begin{array}[]{ll}\\{-\lambda\\}&\text{if
}s\in(0,1)\\\ (-\infty,-\lambda]&\text{if }s=0\\\
{[}-\lambda,+\infty)&\text{if }s=1.\end{array}\right.$

So (6) can be written

$\displaystyle\Phi(x)$ $\displaystyle=\lambda$ $\displaystyle\text{if
}\frac{\mathrm{d}\gamma_{1}}{\mathrm{d}\mu}(x)\in(0,1)$ $\displaystyle\Phi(x)$
$\displaystyle\in[\lambda,+\infty)$ $\displaystyle\text{if
}\frac{\mathrm{d}\gamma_{1}}{\mathrm{d}\mu}(x)=0$ $\displaystyle\Phi(x)$
$\displaystyle\in(-\infty,\lambda]$ $\displaystyle\text{if
}\frac{\mathrm{d}\gamma_{1}}{\mathrm{d}\mu}(x)=1.$

Similarly for $\Psi$. In the discrete case the dual problem and optimality
conditions reduce to the form stated in the proposition. ∎

###  4.1 An Quadratic Time Algorithm

Our algorithm finds the solutions $\gamma$, $\Phi$ and $\Psi$ to the primal-
dual optimality conditions given in Proposition 4.3. We assume that, at
iteration $k$, we have solved
$\mathrm{OPT}(\\{x_{i}\\}_{i=1}^{k-1},\\{y_{j}\\}_{j=1}^{m})$ in the previous
iteration (stored in $L[1:k-1]$) and we now proceed to solve
$\mathrm{OPT}(\\{x_{i}\\}_{i=1}^{k},\\{y_{j}\\}_{j=1}^{m})$ in the current
iteration. Recall that we assume that $\\{x_{i}\\}_{i=1}^{n}$ and
$\\{y_{j}\\}_{j=1}^{m}$ are sorted and that $L[i]$ is the index determining
the transport of mass from $x_{i}$, i.e., if $L[i]\neq-1$ then $x_{i}$ is
associated to $y_{L[i]}$. For simplicity we assume for now that all points in
$\\{x_{i}\\}_{i=1}^{n}$ are distinct, and likewise for
$\\{y_{j}\\}_{j=1}^{m}$. Duplicate points can be handled properly with some
minor additional steps (see appendix). Let $j^{*}$ be the index of the most
attractive $\\{y_{j}\\}_{j=1}^{m}$ for the new point $x_{k}$ under
consideration of the dual variables, i.e.,
$j^{*}=\mathrm{argmin}_{j\in[1:m]}c(x_{k},y_{j})-\Psi_{j}$ and set
$\Phi_{k}=\min\\{c(x_{k},y_{j^{*}})-\Psi_{j^{*}},\lambda\\}$ so that
$\Phi_{k}$ is the largest possible value satisfying the dual constraints, but
no greater than $\lambda$ since at this point the dual objective becomes flat.
We now distinguish three cases:

Case 1:

If $\Phi_{k}=\lambda$, then destroying $x_{k}$ is the most efficient action,
i.e., we set $L[k]=-1$ and proceed to $k+1$.

Case 2:

If $\Phi_{k}<\lambda$ and $y_{j^{*}}$ is unassigned, then we set $L[k]=j^{*}$
and we can proceed to $k+1$.

Case 3:

If $\Phi_{k}<\lambda$ and $y_{j^{*}}$ is already assigned, we must resolve the
conflict between $x_{k}$ and the element currently assigned to $y_{j*}$. This
will be done by a sub-algorithm.

It is easy to see that in the first two cases, if $L$ (or $\gamma$), $\Phi$
and $\Psi$ satisfy the primal-dual conditions of Proposition 4.3 up until
$k-1$, they will now satisfy them up until $k$. Let us now study the third
case in more detail.

One finds that if $y_{j^{*}}$ is already assigned, then it must be to
$x_{k-1}$ (proof in appendix). We now increase $\Phi_{i}$ for $i\in[k-1:k]$
while decreasing $\Psi_{j^{*}}$ (at the same rate) until one of the following
cases occurs:

Case 3.1:

Either $\Phi_{k-1}$ or $\Phi_{k}$ reaches $\lambda$. In this case, the
corresponding $x$ becomes unassigned and the other becomes (or remains)
assigned to $j^{*}$. The conflict is resolved and we proceed to solving the
problem for $k+1$.

Case 3.2:

One reaches the point where $\Phi_{k}+\Psi_{j^{*}+1}=c(x_{k},y_{j^{*}+1})$. In
this case, $x_{k}$ becomes assigned to $y_{j^{*}+1}$, the conflict is resolved
and we move to $k+1$.

Case 3.3a:

One reaches the point where
$\Phi_{k-1}+\Psi_{j^{*}-1}=c(x_{k-1},y_{j^{*}-1})$. If $y_{j^{*}-1}$ is
unassigned, we assign $x_{k-1}$ to $y_{j^{*}-1}$, $x_{k}$ to $y_{j^{*}}$. The
conflict is resolved and we move on.

Case 3.3b:

In the remaining case where $y_{j^{*}-1}$ is already assigned, we will show
that it must be to $x_{k-2}$. This means that the set of points involved in
the conflict increases and we must perform a slight generalization of the
above case 3 iteration until eventually one of the other cases occurs.

At each iteration we consider a contiguous set of
$\\{x_{i}\\}_{i=i_{\min}}^{k-1}$ assigned monotonously to contiguous
$\\{y_{j}\\}_{j=j_{\min}}^{j^{*}}$, where $i_{\min}$ is initially set
$i_{\min}=k-1$ and $j_{\min}$ is the index of $y_{j}$ that has been assigned
to $x_{i_{\min}}$ (i.e., $j_{\min}=j^{*}$), and the additional point $x_{k}$.
We increase $\Phi_{i}$ for $i\in[i_{\min}:k]$ and decrease $\Psi_{j}$ for
$j\in[j_{\min}:j^{\ast}]$, until either $\Phi_{k^{\prime}}$ becomes equal to
$\lambda$ for some $k^{\prime}\in[i_{\min}:k+1]$ (Case 3.1),
$\Phi_{k}+\Psi_{j^{*}+1}=c(x_{k},y_{j^{*}+1})$ (Case 3.2) or
$\Phi_{i_{\min}}+\Psi_{j_{\min}-1}=c(x_{i_{\min}},y_{j_{\min}-1})$ (Case 3.3).
In Cases 3.1 and 3.2 the conflict can be resolved in an obvious way. The same
holds for Case 3.3, when $y_{j_{\min}-1}$ is unassigned (Case 3.3a). Otherwise
(Case 3.3b), one adds one more assigned pair to the chain and restarts the
loop with the pair $(x_{i_{\min}-1},y_{j_{\min}-1})$ added to the chain. Of
course, trivial adaptations have to be made to account for boundary effects,
e.g., Case 3.2 cannot occur if $j^{*}=m$ etcetera. A pseudocode description of
the method is given by Algorithms 1 and 2, and a visual illustration of these
algorithms is provided in Figure 1. Note that for the sake of legibility we
make some simplifications, e.g., boundary checks are ignored (see above), and
we do not specify how to keep track of whether $j^{*}$ is assigned or not.
Also, updating the dual variables at each iteration of the sub-routine yields
a cubic worst case time complexity, but is easier to understand. A complete
version of the algorithm with all checks, appropriate data structures, and
quadratic complexity is given in the appendix. There we also prove the
following claim:

######  Theorem 4.4.

Algorithm 1 is correct, i.e., it is well-defined and returns optimal primal
and dual solutions $L$ and $(\Phi,\Psi)$ to the 1-dimensional OPT problem for
sorted $\\{x_{i}\\}_{i=1}^{n}$, $\\{y_{j}\\}_{j=1}^{m}$. A slight adaptation
of the algorithm (given in the appendix) has a worst case time complexity of
$O(n\max\\{m,n\\})$.

The algorithm can be adjusted to the case where $\lambda_{1}=\infty$ by
setting all instances of $\lambda$ in Algorithms 1 and 2 to $+\infty$ except
for the initialization of $\Psi$. This means that cases 1 and 3.1 never occur.
The algorithm then reduces to a primal-dual version of that given in [4]. If
$\lambda_{2}=\infty$, then one simply flips the two marginals and proceeds as
for $\lambda_{1}=\infty$.

Input: $\\{x_{i}\\}_{i=1}^{n},\\{y_{j}\\}_{j=1}^{m},\lambda$

Output: $L$, $\Psi$, $\Phi$

1 Initialize $\Phi_{i}\leftarrow-\infty$ for $i\in[1:n]$,
$\Psi_{j}\leftarrow\lambda$ for $j\in[1:m]$ and $L_{i}\leftarrow-1$ for
$i\in[1:n]$

2 for  _$k=1,2,\ldots n$_ do

3 $j^{*}\leftarrow\operatorname{argmin}_{j\in[1:m]}c(x_{k},y_{j})-\Psi_{j}$

4 $\Phi_{k}\leftarrow\min\\{c(x_{k},y_{j^{*}})-\Psi_{j^{*}},\lambda\\}$

5 if  _$\Phi_{k}=\lambda$_ then

6 [Case 1] No update on $L$

7 else if  _$j^{*}$ unassigned_ then

8 [Case 2] $L_{k}\leftarrow j^{*}$

9 else

10 [Case 3] Run Algorithm 2.

11

Algorithm 1 opt-1d

Input: ($\\{x_{i}\\}_{i=1}^{n},\\{y_{j}\\}_{j=1}^{m}$, $k$, $j^{*}$, $L$,
$\Phi$, $\Psi$)

Output: (Updated $L$, $\Phi$, $\Psi$, optimal for
$\mathrm{OPT}(\\{x_{i}\\}_{i=1}^{k},\\{y_{j}\\}_{j=1}^{m})$)

1 Initialize $i_{\min}\leftarrow k-1$, $j_{\min}\leftarrow j^{*}$.

2 while  _true_ do

3
$i_{\Delta}\leftarrow\operatorname{argmin}_{i\in[i_{\min}:k]}(\lambda-\Phi_{i})$

4 $\lambda_{\Delta}\leftarrow\lambda-\Phi_{i_{\Delta}}$

5 $\alpha\leftarrow c(x_{k},y_{j^{*}+1})-\Phi_{k}-\Psi_{j^{*}+1}$

6 $\beta\leftarrow
c(x_{i_{\min}},y_{j_{\min}-1})-\Phi_{i_{\min}}-\Psi_{j_{\min}-1}$

7 if  _$\lambda_{\Delta}\leq\min\\{\alpha,\beta\\}$_ then

8 [Case 3.1]

9 $\Phi_{i}\leftarrow\Phi_{i}+\lambda_{\Delta}$ for $i\in[i_{\min}:k]$

10 $\Psi_{j}\leftarrow\Psi_{j}-\lambda_{\Delta}$ for $j\in[j_{\min}:j^{*}]$

11 $L_{i_{\Delta}}\leftarrow-1$, $L_{k}\leftarrow j^{*}$

12 for  _$i\in[i_{\Delta}+1:k-1]$_ do

13 $L_{i}\leftarrow L_{i}-1$

14 return

15

16 else if  _$\alpha\leq\min\\{\lambda_{\Delta},\beta\\}$_ then

17 [Case 3.2]

18 $\Phi_{i}\leftarrow\Phi_{i}+\alpha$ for $i\in[i_{\min}:k]$

19 $\Psi_{j}\leftarrow\Psi_{j}-\alpha$ for $j\in[j_{\min}:j^{*}]$

20 $L_{k}\leftarrow j^{*}+1$

21 return

22 else

23 $\Phi_{i}\leftarrow\Phi_{i}+\beta$ for $i\in[i_{\min}:k]$

24 $\Psi_{j}\leftarrow\Psi_{j}-\beta$ for $j\in[j_{\min}:j^{*}]$

25 if  _$j_{min}-1$ unassigned_ then

26 [Case 3.3a]

27 $L_{i_{\min}}\leftarrow j_{\min}-1$, $L_{k}\leftarrow j^{*}$

28 for  _$i\in[i_{\min}+1:k-1]$_ do

29 $L_{i}\leftarrow L_{i}-1$

30 return

31

32 else

33 [Case 3.3b]

34 $i_{\min}\leftarrow i_{\min}-1$, $j_{\min}\leftarrow j_{\min}-1$

35

36

Algorithm 2 sub-opt Figure 1: Depiction of Algorithms 1 and 2 for solving the
optimal partial transport in one dimension.

###  4.2 Runtime Analysis

We test the wall clock time by sampling the point lists from uniform
distributions and Gaussian mixtures. In particular, we set
$\\{x_{i}\\}_{i=1}^{n}\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\text{
Unif}[-20,20]$,
$\\{y_{j}\\}_{j=1}^{m}\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\text{
Unif}[-40,40]$, $\lambda\in\\{20,100\\}$ and
$\\{x_{i}\\}_{i=1}^{n}\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\frac{1}{5}\sum_{k=1}^{5}\mathcal{N}(-4+2k,1),\\{y_{j}\\}_{j=1}^{m}\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\frac{1}{6}\sum_{k=1}^{6}\mathcal{N}(-5+2k,1)$,
$\lambda\in\\{2,10\\}$, $n\in\\{2000,2500,\ldots 10000\\}$, and $m=n+1000$. We
compare our algorithm with POT (Algorithm 1 in [4]), an unbalanced Sinkhorn
algorithm [11] (we set the entropic regularization parameter to $\lambda/40$),
and linear programming in python OT [22]. Our algorithm, POT and Sinkhorn
algorithm are accelerated by numba and linear programming is written in C++.
Note that POT [4] and the unbalanced Sinkhorn minimize a different model. In
addition, for the latter the performance depends strongly on the strength of
regularization and we found that it was not competitive in the regime of low
blur. For each $(n,m)$, we repeat the computation 10 times and compute the
average. For our method and POT, the time of sorting is included, and for the
linear programming and Sinkhorn algorithms, the time of computing the cost
matrix is included. We also visualize our algorithm’s solutions for
$\\{x_{i}\\}_{i=1}^{n}\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\text{Unif}(-20,20),\\{y_{j}\\}_{j=1}^{m}\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\text{Unif}(-40,40),n=8,m=16$
and $\lambda\in\\{1,10,100,1000\\}$ (see figure 3). One can observe that as
$\lambda$ increases larger transportations are permitted. The data type is
32-bit float number and the experiments are conducted on a Linux computer with
AMD EPYC 7702P CPU with 64 cores and 256GB DDR4 RAM.

(a) wall-clock time for uniform distributions

(b) Wall-clock time for Gaussian mixture distribution

Figure 2: We test the wall-clock time of our Algorithm 1, the partial OT
solver (Algorithm 1 in [4]), the unbalanced Sinkhorn algorithm [11], and the
linear programming solver in POT [22], which is written in C. The maximum
number of iterations for linear programming is set to $100n$ and it has been
reached several times in the experiment. Figure 3: Output of the proposed
algorithm for $\lambda\in[1,10,100,1000]$ on a sample pair of measures.

##  5 Sliced Optimal Partial Transport

In practice, data has multiple dimensions and the 1-D computational methods
cannot be applied directly. In the balanced OT setting the sliced OT approach
[45, 5, 31, 33, 30] applies the 1-D OT solver on projections (i.e., slices) of
two r-dimensional distributions. Inspired by these works, we extend the sliced
OT technique into the OPT setting and introduce the so called “sliced-
unbalanced optimal transport” problem.

######  Definition 5.1.

In $\mathbb{R}^{r}$ space, given $\mu,\nu\in\mathcal{M}_{+}(\Omega)$ where
$\Omega\subset\mathbb{R}^{r}$ and $\lambda:\mathbb{S}^{r-1}\to\mathbb{R}_{++}$
is an $\mathrm{L}^{1}$ function, we define the sliced optimal partial
transport (SOPT) problem as follows:

$\displaystyle\mathrm{SOPT}_{\lambda}(\mu,\nu)=\int_{\mathbb{S}^{r-1}}\mathrm{OPT}_{\lambda(\theta)}(\langle\theta,\cdot\rangle_{\\#}\mu,\langle\theta,\cdot\rangle_{\\#}\nu)\,\mathrm{d}\sigma(\theta)$
(7)

where $\mathrm{OPT}_{\lambda}(\cdot,\cdot)$ is defined in (3),
$\sigma\in\mathcal{P}(\mathbb{S}^{r-1})$ is a probability measure such that
$\mathrm{supp}(\sigma)=\mathbb{S}^{r-1}$.

Generally $\sigma$ is set to be the uniform distribution on the unit ball
$\mathbb{S}^{r-1}$. Whenever $\mathrm{supp}(\sigma)=\mathbb{S}^{r-1}$,
$\mathrm{SOPT}_{\lambda}(\mu,\nu)$ defines a metric.

######  Theorem 5.2.

Suppose $c:\Omega\times\Omega\to\mathbb{R}_{+}$ is a metric on $\Omega$ and
$\lambda\in\mathrm{L}^{1}(\mathbb{S}^{r-1};\mathbb{R}_{++})$ satisfies
$\inf_{\theta\in\mathbb{S}^{r-1}}\lambda(\theta)>0$, then
$\mathrm{SOPT}_{\lambda}(\mu,\nu)$ is a metric in $\mathcal{M}_{+}(\Omega)$.

In practice, this integration is usually approximated using a Monte Carlo
scheme that draws a finite number of i.i.d. samples
$\\{\theta_{l}\\}_{l=1}^{N}$ from $\text{Unif}(\mathbb{S}^{d-1})$ and replaces
the integral with an emprical average:

$\mathrm{SOPT}_{\lambda}(\mu,\nu)\approx\frac{1}{N}\sum_{l=1}^{N}\mathrm{OPT}_{\lambda_{l}}(\langle\theta_{l},\cdot\rangle_{\\#}\mu,\langle\theta_{l},\cdot\rangle_{\\#}\nu).$

###### Proof of Theorem 5.2.

First we claim
$\mathrm{SOPT}_{\lambda}(\cdot,\cdot):(\mathcal{M}_{+}(\Omega))^{2}\to\mathbb{R}_{+}$
is a well defined function. It is clear $\mathrm{SOPT}_{\lambda}(\cdot,\cdot)$
is a function with domain $\mathcal{M}_{+}(\Omega)^{2}$ and co-domain
$\mathbb{R}\cup\\{\pm\infty\\}$. Pick $\mu,\nu$, we will show
$\mathrm{SOPT}_{\lambda}\in[0,\infty)$. We have

$\displaystyle\mathrm{SOPT}_{\lambda}(\mu,\nu)$
$\displaystyle=\int_{\mathbb{S}^{d-1}}\mathrm{OPT}_{\lambda(\theta)}(\langle\theta,\cdot\rangle_{\\#}\mu,(\langle\theta,\cdot\rangle_{\\#}\nu)\,\mathrm{d}\sigma(\theta)\geq
0$ (8)

where the inequality follows from the fact
$\mathrm{OPT}{\lambda(\theta)}(\langle\theta,\cdot\rangle_{\\#}\mu,(\langle\theta,\rangle_{\\#}\nu)\geq
0,\forall\theta$. It remains to show
$\mathrm{SOPT}_{\lambda}(\mu,\nu)<\infty$. We have

$\displaystyle\mathrm{SOPT}_{\lambda}(\mu,\nu)$
$\displaystyle\leq\int_{\mathbb{S}^{d-1}}\frac{1}{2}\lambda(\theta)(\|\langle\theta,\cdot\rangle_{\\#}\mu\|_{\mathrm{TV}}+\|\langle\theta,\cdot\rangle_{\\#}\nu\|_{\mathrm{TV}})\,\mathrm{d}\sigma(\theta)$
$\displaystyle=\frac{1}{2}(\mu(\Omega)+\nu(\Omega))\int_{\mathbb{S}^{d-1}}\lambda(\theta)\,\mathrm{d}\sigma(\theta)$
$\displaystyle<\infty$

where the first inequality follows by plugging zero measure into the cost
function in (3); the second inequality holds since $\lambda$ is an
$\mathrm{L}^{1}$ function.

Next, we will show $\mu=\nu$ iff $\mathrm{SOPT}_{\lambda}(\mu,\nu)=0$. If
$\mu=\nu$, we have for every $\theta$,
$\langle\theta,\rangle_{\\#}\mu=\langle\theta,\rangle_{\\#}\nu$ and thus
$\mathrm{OPT}_{\lambda(\theta)}(\langle\theta,\cdot\rangle_{\\#}\mu,\langle\theta,\cdot\rangle_{\\#}\nu)=0$.
Therefore $\mathrm{SOPT}_{\lambda}(\mu,\nu)=0$. For the reverse direction, we
suppose $\mathrm{SOPT}_{\lambda}(\mu,\nu)=0$. Since
$\mathrm{supp}(\sigma)=\mathbb{S}^{r-1}$, we have for almost every $\theta$,
$\mathrm{OPT}_{\lambda(\theta)}(\langle\theta,\cdot\rangle_{\\#}\mu,\langle\theta,\cdot\rangle_{\\#}\nu)=0$.
For almost every $\theta$, since $\lambda(\theta)>0$, and
$\mathrm{OPT}_{\lambda(\theta)}(\cdot,\cdot)$ is a metric (see [12,
Proposition 2.10] or [43, Proposition 5]), we have
$\langle\theta,\cdot\rangle_{\\#}\mu=\langle\theta,\cdot\rangle_{\\#}\nu$. By
the injectivity of Radon transform on measures (see [5, Proposition 7]), we
have $\mu=\nu$.

For symmetry, we have

$\displaystyle\mathrm{SOPT}_{\lambda}(\mu,\nu)$
$\displaystyle=\int_{\mathbb{S}^{d-1}}\mathrm{OPT}_{\lambda(\theta)}(\langle\theta,\cdot\rangle_{\\#}\mu,\langle\theta,\rangle_{\\#}\nu)\,\mathrm{d}\sigma(\theta)$
$\displaystyle=\int_{\mathbb{S}^{d-1}}\mathrm{OPT}_{\lambda(\theta)}(\langle\theta,\cdot\rangle_{\\#}\nu,\langle\theta,\cdot\rangle_{\\#}\mu)\,\mathrm{d}\sigma(\theta)$
$\displaystyle=\mathrm{SOPT}_{\lambda}(\nu,\mu)$ (9)

where the second equality follows from the fact
$\mathrm{OPT}_{\lambda(\theta)}(\cdot,\cdot)$ is a metric for each $\theta$.

For the triangle inequality, we choose $\mu_{1},\mu_{2},\mu_{3}$ and have:

$\displaystyle\mathrm{SOPT}_{\lambda}(\mu_{1},\mu_{3})$
$\displaystyle=\int_{\mathbb{S}^{d-1}}\mathrm{OPT}_{\lambda(\theta)}(\langle\theta,\cdot\rangle_{\\#}\mu_{1},\langle\theta,\cdot\rangle_{\\#}\mu_{3})\,\mathrm{d}\sigma(\theta)$
$\displaystyle\leq\int_{\mathbb{S}^{d-1}}\left[\mathrm{OPT}_{\lambda(\theta)}(\langle\theta,\cdot\rangle_{\\#}\mu_{1},\langle\theta,\cdot\rangle_{\\#}\mu_{2})+\mathrm{OPT}_{\lambda(\theta)}(\langle\theta,\cdot\rangle_{\\#}\mu_{2},\langle\theta,\cdot\rangle_{\\#}\mu_{3})\right]\,\mathrm{d}\sigma(\theta)$
$\displaystyle=\left(\int_{\mathbb{S}^{d-1}}\mathrm{OPT}_{\lambda(\theta)}(\langle\theta,\cdot\rangle_{\\#}\mu_{1},\langle\theta,\cdot\rangle_{\\#}\mu_{2})\,\mathrm{d}\sigma(\theta)\right)+\left(\int_{\mathbb{S}^{d-1}}\mathrm{OPT}_{\lambda(\theta)}(\langle\theta,\cdot\rangle_{\\#}\mu_{1},\langle\theta,\cdot\rangle_{\\#}\mu_{2})\,\mathrm{d}\sigma(\theta)\right)$
$\displaystyle=\mathrm{SOPT}_{\lambda}(\mu_{1},\mu_{2})+\mathrm{SOPT}_{\lambda}(\mu_{2},\mu_{3})$

where the first inequality follows from the fact
$\mathrm{OPT}_{\lambda(\theta)}(\cdot,\cdot)$ is a metric for each $\theta$. ∎

##  6 Applications

###  6.1 Point Cloud Registration

Point cloud registration is a transformation estimation problem between two
point clouds, which is a critical problem in numerous computer vision
applications [51, 27]. In particular, given two 3-D point clouds,
$x_{i}\sim\mu$, $i=1,\dots,n$ and $y_{j}\sim\nu$, $j=1,\dots,m$, one assumes
there is an unknown mapping, $T$, that satisfies $\nu=T_{\\#}\mu$. In many
applications, the mapping $T$ is restricted to have the following form,
$Tx=sRx+\beta$, where $R$ is a $3\times 3$ dimensional rotation matrix, $s>0$
is the scaling and $\beta$, called translation vector, is $3\times 1$ vector.
The goal is then to estimate the transform, $T$, based on the two point
clouds.

| 10k,5% | 10k,7% | 9k,5% | 9k,7%
---|---|---|---|---
ICP-D | 1.10(1.59) | 3.60(0.11) | 1.65(1.13) | 2.04(2.28)
ICP-U | 3.72(0.53) | 3.72(0.52) | 3.69(0.63) | 3.92(0.32)
SPOT | 1.27(0.01) | 1.40(0.15) | 1.42(0.13) | 1.53(1.50)
Ours | 0.01(1e-3) | 0.02(2e-3) | 0.20(0.09) | 0.33(0.03)
Table 1: We compute the mean (and variance, in parenthesis) of errors in the Frobenius norm between the ground truth and estimated transportation matrices for ICP(Du), ICP(Umeyama), SPOT and our method. We vary the size of the source point cloud from 9k to 10k samples, the percentage of noise from 5% to 7% (on both source and target datasets). Figure 4: We visualize the results of ICP(Du), SPOT and our method. In each image, the target point cloud is in green and source point cloud is in orange. |  |  | ICP | SPOT | Ours
---|---|---|---|---|---
Dataset | Source | Target | Time/iter. | # Iter. | Time/iter. | # Iter. | Time/iter. | # Iter.
Bunny (5%) | $9k$ | $10k$ | $0.66s$ | $50$ | $10.9s$ | $100$ | $0.31s$ | $2000$
Bunny (7%) | $10k$ | $10k$ | $0.76s$ | $130$ | $13.2s$ | $100$ | $0.35s$ | $2000$
Dragon (5%) | $9k$ | $10k$ | $0.66s$ | $100$ | $10.4s$ | $100$ | $0.31s$ | $2000$
Dragon (7%) | $10k$ | $10k$ | $0.77s$ | $100$ | $13.1s$ | $80$ | $0.35s$ | $1500$
Mumble (5%) | $9k$ | $10k$ | $0.65s$ | $100$ | $10.78s$ | $100$ | $0.32s$ | $2000$
Mumble (7%) | $10k$ | $10k$ | $0.78s$ | $100$ | $13.18s$ | $80$ | $0.36s$ | $1500$
Castle (5%) | $9k$ | $10k$ | $0.66s$ | $150$ | $10.7s$ | $100$ | $0.31s$ | $2000$
Castle (7%) | $10k$ | $10k$ | $0.76s$ | $350$ | $13.7s$ | $80$ | $0.35s$ | $1800$
Table 2: This table reports the data for our method in shape registration
experiment: the number of source and target distribution, the percentage of
noise, the computation times per iteration and numbers of iterations they took
to converge for ICP(Du), SPOT, and our method. The source point cloud is in
orange color and the target is in blue. The percentage of noise is given in
brackets next to the dataset name.

The classic approach for solving this problem is Iterative Closest Point
Algorithms (ICP) introduced by [9, 3]. By [53]’s work, classical ICP can be
extended into the uniformly scaled setting, with further developments by [17].
To address the some issues of ICP methods (convergence to local minimum, poor
performance when the size of the two datasets are not equal), [4] proposed the
Fast Iterative Sliced Transport algorithm (FIST) using sliced partial optimal
transport (SPOT).

Problem setup and our method. We consider the uniform scaled point cloud
registration problem and assume a subset of points in both the source and
target point clouds are corrupted with additional uniformly distributed noise.
We suppose prior knowledge of the proportion of noise is given (i.e., we have
prior knowledge of number of clean data).

In general, the registration problem can be iteratively solved and each
iteration contains two steps: estimating the correspondence and computing the
optimal transform from corresponding points. The second step has a closed form
solution. For the first step, the ICP method estimates the correspondence by
finding the closest $y$ for each transformed $x$. Inspired by this work, we
estimate the correspondence by using our SOPT solver. See Algorithm 3.

Input: $\\{x_{i}\\}_{i=1}^{n},\\{y_{j}\\}_{j=1}^{m},$ $n_{0}$:the # of clean
$x$ , $N$: # of projections

Output: $R,s,\beta$

1 initialize $R,s,\beta,\lambda$, sample
$\\{\theta_{i}\\}_{i=1}^{N}\subset\mathbb{S}^{2}$

2 for  _$l=1,\ldots N$_ do

3 $\hat{Y}\leftarrow sRX+\beta$

4 Compute transportation plan $L$ of
$\mathrm{OPT}_{\lambda}(\theta_{l}^{T}\hat{Y},\theta_{l}^{T}Y)$ by algorithm 1

5 $\forall i\in\mathrm{dom}(L)$,
$\hat{y}_{i}\leftarrow\hat{y}_{i}+(\theta_{l}^{T}y_{L[i]}-\theta_{l}^{T}\hat{y}_{i})\theta$

6 Compute $R,s,\beta$ from $(X[dom(L)],\hat{Y}[dom(L)])$ by ICP (e.g.,
equations (39)-(42) in [53])

7 If $|dom(L)|>n_{0}$, decrease $\lambda$; otherwise, increase $\lambda$.

8

Algorithm 3 iterative-sopt

Experiment. We illustrate our algorithm on different 3D point clouds,
including Stanford Bunny, Stanford dragon, Witch-castle and Mumble Sitting.
For each dataset, we generate transforms by uniformly sampling angles from
$[-1/3\pi,1/3\pi]$, translations from $[-2std,2std]$, and scalings from
$[0,2]$, where $std$ is the standard deviation of the dataset. Then we sample
noise uniformly from the region $[-M,M]^{3}$ where
$M=\max_{i\in[1:n]}(\|x_{i}\|)$ and concatenate it to our point clouds. The
number of points in the target (clean) data is fixed to be $10k$ and we vary
the number of points in the source (clean) data from $9k$ to $10k$, and the
percentage of noise from $5\%$ to $7\%$.

Performance. For the accuracy, we compute the average and standard deviation
of error defined by the Frobenius norm between the estimated transformation
and the ground truth (see Table 1). For the accuracy, we observe ICP methods
systematically fails, since the nearest neighbor matching technique induces an
non-injective correspondence, which may result in a too small scaling factor.
SPOT can successfully recover the rotation, but it fails to recover the
scaling. Our method is the only one that recovers the ground truth for the
noise corrupted data (since it utilizes the prior knowledge).

For running time, ICP has the fastest convergence time which is generally
100-260 second since finding the correspondence by closest neighbor can be
done immediately after the cost matrix is computed. SPOT requires 1000-1300
second and our method requires 500-700 seconds. The data type is $32$-bit
float number and the experiment is conducted on a Linux computer with AMD EPYC
7702P CPU with 64 cores and 256GB DDR4 RAM.

Figure 5: We visualize the processed point cloud for every methods with
respect to time. The dataset is Stanford dragon. In each image, the point
cloud in orange is source and the point cloud in green color is target.

###  6.2 Color Adaptation

Transferring colors between images is a classical task in computer vision and
image science. Given two images, the goal is to impose on one of the images
(source image) the histogram of the other image (target image). Optimal
transport-based approaches have been developed and achieved great success in
this task [11, 5, 45, 19]. However, in the balanced OT setting, the OT-based
approach requires normalizing the histograms of (sampled) colors, which may
lead to undesired performance. For example, suppose the majority of a target
image is red, (e.g. an image of evening sunset) and the majority of a source
image is green (e.g. image of trees). Then balanced-OT-based approaches will
produce a red tree in the result. To address this issue, [4] applied the SPOT-
based approach which will match all the pixels in the source image to partial
pixels in the target image.

Our method. Inspired by [19, 4], our method contains the following three
steps: First, we sample pixels from the source and target image by k-mean
clustering (or another sampling method). Second, we transport the sampled
source pixels into the target domain. If OT or entropic OT is applied, it can
be done by the optimal transportation plan; if sliced-OT is applied, the
source pixels would be updated iteratively for each slice based on the
gradient of 1-D OT with respect to the source pixels (see equation (46) in
[5], or line 5 in our algorithm 3). In our method, we apply the transportation
plan from OPT. Third, we reconstruct the source image based on the transported
source pixels (e.g. see Equation 4.1 in [19]).

Experiment. We first normalize all the pixels in the source and target images
to be in range $[0,1]$, then we use k-means clustering to sample $5000$ pixels
from the source image and $10000$ pixels from the target image. We compare the
performance of the OT-based and Entropic-OT based domain adaptation functions
in PythonOT [22] (ot.da.EMDTransport and ot.da.SinkhornTransport) whose OT’s
solver is written in C++ ***We modify their code to increase the speed., SPOT
[4] and our method based on sliced optimal partial transport. For our method,
we test it in two schemes, $\lambda=10.0$ and $\lambda<2.0$. In the first
case, $\lambda$ achieves the maximum of distance of two (normalized) pixels,
that is, we will transport all the source pixels into target domain. In second
case, we choose $\lambda$ and theoretically, only partial source pixels will
be transported into the target domain.

Performance. In these examples, the OT-based approach which matches all
(sampled) pixels of the source image to all pixels of a target image can lead
to undesired results. For example, in the second row of Figure 6, the third
image has dark blue on the sky and red color on the ground. This issue is
alleviated in SPOT and our method. In our method, when $\lambda=5.0$, we will
transfer all the (sampled) pixels from source to target and the result is
similar to the result of SPOT ***We conjecture the two results are not exactly
the same due to the randomness of projections.. When $\lambda<2.0$, the result
image is closer to the source image. OT-based method requires 40-50 seconds
(we set the maximum iteration number of linear programming to be 1000,000);
Partial OT method requires 80-90 seconds (the # of projections is set to be
400) and our method requires 60-80 seconds (the # of projections is set to be
400). The data type is 32-bit float number and the experiment is conducted on
a Linux computer with AMD EPYC 7702P CPU with 64 cores and 256GB DDR4 RAM.

Figure 6: We transfer colors from source image to the target image by the
methods based on optimal transport [19], SPOT [4] and our SOPT. For our
method, we set $\lambda=10$ and a small value less than $2$. Image via Flickr:
Facade by Phil Whitehouse, palace by Neil Williamson, clouds by Tim Wang, air
balloon by Kirt Edblom, roses by Felix Schaumburg.

##  7 Conclusion and future work

This paper proposes a fast computational method for solving the OPT problem
for one-dimensional discrete measures. We provide computational and wall-clock
analysis experiments to assess our proposed algorithm’s correctness and
computational benefits. Then, utilizing one-dimensional slices of an
$r$-dimensional measure we propose the “sliced optimal partial transport
(SOPT)” distance. Beyond our theoretical and algorithmic contributions, we
provide empirical evidence that SOPT is practical for large-scale applications
like point cloud registration and image color adaptation. In point cloud
registration, we show that compared with other classical methods, our SOPT-
based approach adds robustness when the source and target point clouds are
corrupted by noise. In the future, we will investigate the potential
applications of SOPT in other machine learning tasks such as open set domain
adaptation problems and measuring task similarities in continual and
curriculum learning.

##  8 Acknowledgment

The authors thank Rana Muhammad Shahroz Khan for helping in code testing and
thank Dr. Hengrong Du<EMAIL_ADDRESS>for helpful discussions.
This work was partially supported by the Defense Advanced Research Projects
Agency (DARPA) under Contract No. HR00112190132. MT was supported by the
European Research Council under the European Union’s Horizon 2020 research and
innovation programme Grant Agreement No. 777826 (NoMADS).

## References

  * [1] Martin Arjovsky, Soumith Chintala, and Léon Bottou.  Wasserstein generative adversarial networks.  In International conference on machine learning, pages 214–223. PMLR, 2017.
  * [2] Jean-David Benamou, Guillaume Carlier, Marco Cuturi, Luca Nenna, and Gabriel Peyré.  Iterative bregman projections for regularized transportation problems.  SIAM Journal on Scientific Computing, 37(2):A1111–A1138, 2015.
  * [3] Paul J Besl and Neil D McKay.  Method for registration of 3-d shapes.  In Sensor fusion IV: control paradigms and data structures, volume 1611, pages 586–606. Spie, 1992.
  * [4] Nicolas Bonneel and David Coeurjolly.  SPOT: sliced partial optimal transport.  ACM Transactions on Graphics, 38(4):1–13, 2019.
  * [5] Nicolas Bonneel, Julien Rabin, Gabriel Peyré, and Hanspeter Pfister.  Sliced and Radon Wasserstein barycenters of measures.  Journal of Mathematical Imaging and Vision, 51(1):22–45, 2015.
  * [6] R. E. Burkhard, B. Klinz, and R. Rudolf.  Perspectives of Monge properties in optimization.  Discr. Appl. Math., 70(2):95–161, 1996.
  * [7] Luis A Caffarelli and Robert J McCann.  Free boundaries in optimal transport and monge-ampere obstacle problems.  Annals of mathematics, pages 673–730, 2010.
  * [8] Yongxin Chen, Tryphon T Georgiou, Lipeng Ning, and Allen Tannenbaum.  Matricial wasserstein-1 distance.  IEEE control systems letters, 1(1):14–19, 2017.
  * [9] Yang Chen and Gérard Medioni.  Object modelling by registration of multiple range images.  Image and vision computing, 10(3):145–155, 1992.
  * [10] Lenaic Chizat, Gabriel Peyré, Bernhard Schmitzer, and François-Xavier Vialard.  An interpolating distance between optimal transport and Fisher–Rao metrics.  Foundations of Computational Mathematics, 18(1):1–44, 2018.
  * [11] 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.
  * [12] Lenaic Chizat, Gabriel Peyré, Bernhard Schmitzer, and François-Xavier Vialard.  Unbalanced optimal transport: Dynamic and Kantorovich formulations.  Journal of Functional Analysis, 274(11):3090–3123, 2018.
  * [13] Nicolas Courty, Rémi Flamary, Amaury Habrard, and Alain Rakotomamonjy.  Joint distribution optimal transportation for domain adaptation.  Advances in Neural Information Processing Systems, 30, 2017.
  * [14] 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, pages 274–289. Springer, 2014.
  * [15] Marco Cuturi.  Sinkhorn distances: Lightspeed computation of optimal transport.  Advances in neural information processing systems, 26, 2013.
  * [16] Ishan Deshpande, Yuan-Ting Hu, Ruoyu Sun, Ayis Pyrros, Nasir Siddiqui, Sanmi Koyejo, Zhizhen Zhao, David Forsyth, and Alexander G Schwing.  Max-sliced wasserstein distance and its use for gans.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 10648–10656, 2019.
  * [17] Shaoyi Du, Nanning Zheng, Shihui Ying, Qubo You, and Yang Wu.  An extension of the icp algorithm considering scale factor.  In 2007 IEEE International Conference on Image Processing, volume 5, pages V–193. IEEE, 2007.
  * [18] Kilian Fatras, Thibault Séjourné, Rémi Flamary, and Nicolas Courty.  Unbalanced minibatch optimal transport; applications to domain adaptation.  In International Conference on Machine Learning, pages 3186–3197. PMLR, 2021.
  * [19] Sira Ferradans, Nicolas Papadakis, Gabriel Peyré, and Jean-François Aujol.  Regularized discrete optimal transport.  SIAM Journal on Imaging Sciences, 7(3):1853–1882, 2014.
  * [20] Alessio Figalli.  The optimal partial transport problem.  Archive for rational mechanics and analysis, 195(2):533–560, 2010\.
  * [21] Alessio Figalli and Nicola Gigli.  A new transportation distance between non-negative measures, with applications to gradients flows with dirichlet boundary conditions.  Journal de mathématiques pures et appliquées, 94(2):107–130, 2010.
  * [22] Rémi Flamary, Nicolas Courty, Alexandre Gramfort, Mokhtar Z. Alaya, Aurélie Boisbunon, Stanislas Chambon, Laetitia Chapel, Adrien Corenflos, Kilian Fatras, Nemo Fournier, Léo Gautheron, Nathalie T.H. Gayraud, Hicham Janati, Alain Rakotomamonjy, Ievgen Redko, Antoine Rolet, Antony Schutz, Vivien Seguy, Danica J. Sutherland, Romain Tavenard, Alexander Tong, and Titouan Vayer.  Pot: Python optimal transport.  Journal of Machine Learning Research, 22(78):1–8, 2021.
  * [23] Charlie Frogner, Chiyuan Zhang, Hossein Mobahi, Mauricio Araya, and Tomaso A Poggio.  Learning with a Wasserstein loss.  Advances in neural information processing systems, 28, 2015.
  * [24] Aude Genevay, Marco Cuturi, Gabriel Peyré, and Francis Bach.  Stochastic optimization for large-scale optimal transport.  Advances in neural information processing systems, 29, 2016.
  * [25] Aude Genevay, Gabriel Peyré, and Marco Cuturi.  GAN and VAE from an optimal transport point of view.  arXiv preprint arXiv:1706.01807, 2017.
  * [26] Kevin Guittet.  Extended Kantorovich norms: a tool for optimization.  PhD thesis, INRIA, 2002.
  * [27] Xiaoshui Huang, Guofeng Mei, Jian Zhang, and Rana Abbas.  A comprehensive survey on point cloud registration.  arXiv preprint arXiv:2103.02690, 2021.
  * [28] Leonid V Kantorovich.  On a problem of monge.  In CR (Doklady) Acad. Sci. URSS (NS), volume 3, pages 225–226, 1948\.
  * [29] Narendra Karmarkar.  A new polynomial-time algorithm for linear programming.  In Proceedings of the sixteenth annual ACM symposium on Theory of computing, pages 302–311, 1984.
  * [30] Soheil Kolouri, Kimia Nadjahi, Umut Simsekli, Roland Badeau, and Gustavo Rohde.  Generalized sliced Wasserstein distances.  Advances in Neural Information Processing Systems, 32, 2019.
  * [31] Soheil Kolouri, Se Rim Park, and Gustavo K Rohde.  The radon cumulative distribution transform and its application to image classification.  IEEE transactions on image processing, 25(2):920–934, 2015.
  * [32] Soheil Kolouri, Se Rim Park, Matthew Thorpe, Dejan Slepcev, and Gustavo K Rohde.  Optimal mass transport: Signal processing and machine-learning applications.  IEEE signal processing magazine, 34(4):43–59, 2017.
  * [33] Soheil Kolouri, Yang Zou, and Gustavo K Rohde.  Sliced wasserstein kernels for probability distributions.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 5258–5267, 2016.
  * [34] Tam Le, Makoto Yamada, Kenji Fukumizu, and Marco Cuturi.  Tree-sliced variants of wasserstein distances.  Advances in neural information processing systems, 32, 2019.
  * [35] Wonjun Lee, Rongjie Lai, Wuchen Li, and Stanley Osher.  Generalized unnormalized optimal transport and its fast algorithms.  Journal of Computational Physics, 436:110041, 2021.
  * [36] Jan Lellmann, Dirk A Lorenz, Carola Schonlieb, and Tuomo Valkonen.  Imaging with kantorovich–rubinstein discrepancy.  SIAM Journal on Imaging Sciences, 7(4):2833–2859, 2014.
  * [37] Matthias Liero, Alexander Mielke, and Giuseppe Savare.  Optimal entropy-transport problems and a new Hellinger–Kantorovich distance between positive measures.  Inventiones mathematicae, 211(3):969–1117, 2018.
  * [38] Huidong Liu, Xianfeng Gu, and Dimitris Samaras.  Wasserstein GAN with quadratic transport cost.  In Proceedings of the IEEE/CVF international conference on computer vision, pages 4832–4841, 2019.
  * [39] Antoine Liutkus, Umut Simsekli, Szymon Majewski, Alain Durmus, and Fabian-Robert Stöter.  Sliced-wasserstein flows: Nonparametric generative modeling via optimal transport and diffusions.  In International Conference on Machine Learning, pages 4104–4113. PMLR, 2019.
  * [40] Grégoire Montavon, Klaus-Robert Müller, and Marco Cuturi.  Wasserstein training of restricted Boltzmann machines.  Advances in Neural Information Processing Systems, 29, 2016.
  * [41] James Orlin.  A faster strongly polynomial minimum cost flow algorithm.  In Proceedings of the Twentieth annual ACM symposium on Theory of Computing, pages 377–387, 1988.
  * [42] Gabriel Peyré, Marco Cuturi, et al.  Computational optimal transport: With applications to data science.  Foundations and Trends® in Machine Learning, 11(5-6):355–607, 2019.
  * [43] 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\.
  * [44] Benedetto Piccoli and Francesco Rossi.  On properties of the generalized wasserstein distance.  Archive for Rational Mechanics and Analysis, 222(3):1339–1365, 2016\.
  * [45] Julien Rabin, Gabriel Peyré, Julie Delon, and Marc Bernot.  Wasserstein barycenter and its application to texture mixing.  In International Conference on Scale Space and Variational Methods in Computer Vision, pages 435–446. Springer, 2011.
  * [46] Filippo Santambrogio.  Optimal Transport for Applied Mathematicians, volume 87 of Progress in Nonlinear Differential Equations and Their Applications.  Birkhäuser Boston, 2015.
  * [47] Ryoma Sato, Marco Cuturi, Makoto Yamada, and Hisashi Kashima.  Fast and robust comparison of probability measures in heterogeneous spaces.  arXiv preprint arXiv:2002.01615, 2020.
  * [48] Richard Sinkhorn.  A relationship between arbitrary positive matrices and doubly stochastic matrices.  The annals of mathematical statistics, 35(2):876–879, 1964.
  * [49] 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.
  * [50] Justin Solomon, Raif Rustamov, Leonidas Guibas, and Adrian Butscher.  Wasserstein propagation for semi-supervised learning.  In International Conference on Machine Learning, pages 306–314. PMLR, 2014.
  * [51] Richard Szeliski.  Computer vision: algorithms and applications.  Springer Science & Business Media, 2010.
  * [52] Ilya Tolstikhin, Olivier Bousquet, Sylvain Gelly, and Bernhard Schoelkopf.  Wasserstein auto-encoders.  arXiv preprint arXiv:1711.01558, 2017.
  * [53] Shinji Umeyama.  Least-squares estimation of transformation parameters between two point patterns.  IEEE Transactions on Pattern Analysis & Machine Intelligence, 13(04):376–380, 1991.
  * [54] Cedric Villani.  Optimal transport: old and new.  Springer, 2009.

##  Appendix A Relation Between Optimal Partial Transport and Unbalanced
Optimal Transport

If we replace the penalty term of (3) with a constraint, i.e. we impose the
condition $\gamma(\Omega)\geq M$, and use the fact that
$|\pi_{1\\#}\gamma|=|\pi_{2\\#}\gamma|=\gamma(\Omega)$, then (3) is closely
related to the Lagrangian formulation of the following “primal problem”:

$\displaystyle\text{Primal-OPT}(\mu,\nu;M)$
$\displaystyle=\inf_{\gamma\in\mathcal{M}_{+}(\Omega^{2})}\int
c(x,y)\,\mathrm{d}\gamma(x,y)\quad\text{s.t. }\gamma(\Omega)\geq M.$ (10)

This, in turn, is closely related to the optimal partial transport problem
proposed by [7, 21, 20] (the difference being the mass constraint of $\gamma$
is imposed as an equality $\gamma(\Omega)=M$ rather than a lower bound
$\gamma(\Omega)\geq M$).

Another equivalent form of the OPT problem defined in (3) is the “generalized
Wasserstein distance” in [43, 44] (We refer to [10, Proposition 1.1] and [43,
Proposition 4] for their equivalence.) Recently, more systematic studies of
so-called “unbalanced optimal transport” or “optimal entropy transportation”
problems have been conducted, for instance in [12] and [37]. OPT, (3), can be
seen as a special case of these models, see for instance [12, Theorem 5.2]. It
is also well known that in addition to the static Kantorovich formulations
presented here, one can also give equivalent dynamic formulations in the
spirit of the Benamou–Brenier formula, e.g. [12]. Finally, a related class of
models with close relations to the POT problem is discussed in [35] under the
name “Generalized Unnormalized Optimal Transport”(GUOT).

##  Appendix B Relation Between Optimal Partial Transport and Optimal
Transport

Inspired by Caffarelli et al.s’ technique [7], suppose
$\Omega=\mathbb{R}^{r}$, we introduce an isolated point $\hat{\infty}$ into
$\Omega$ by letting $\hat{\Omega}=\Omega\cup\\{\hat{\infty}\\}$. Suppose
$\hat{\mu}=\mu+\nu(\Omega)\delta_{\hat{\infty}}$, where $\delta_{\hat{x}}$ is
the Dirac mass at $\hat{x}\in\hat{\Omega}$ and
$\hat{\nu}=\nu+\mu(\Omega)\delta_{\hat{\infty}}$,
$\hat{c}(x,y):\hat{\Omega}\times\hat{\Omega}\to\mathbb{R}_{+}$ is defined as

$\hat{c}(x,y):=\begin{cases}c(x,y)-\lambda&\text{if }x,y\neq\hat{\infty}\\\
0&\text{otherwise.}\end{cases}$

We introduce the following optimal transport problem:

$\displaystyle\inf_{\hat{\gamma}\in\Gamma(\hat{\mu},\hat{\nu})}\int\hat{c}(x,y)\,\mathrm{d}\hat{\gamma}(x,y)$
(11)

We claim there exist an equivalence between this OT problem and
$\mathrm{OPT}\lambda(\mu,\nu)$.

######  Proposition B.1.

The mapping: $T:\Gamma_{\leq}(\mu,\nu)\to\Gamma(\hat{\mu},\hat{\nu})$ defined
by

$\gamma\mapsto\hat{\gamma}=\gamma+(\mu-(\pi_{1})_{\\#}\gamma)\otimes\delta_{\hat{\infty}}+\delta_{\hat{\infty}}\otimes(\nu-(\pi_{2})_{\\#}\gamma)+\gamma(\Omega^{2})\delta_{\hat{\infty},\hat{\infty}},$
(12)

is a bijection and $\gamma$ is optimal in $\mathrm{OPT}_{\lambda}(\mu,\nu)$ if
and only if $\hat{\gamma}$ is optimal in (11).

###### Proof.

First we will show that $\hat{\gamma}=T(\gamma)\in\Gamma(\hat{\mu},\hat{\nu})$
for $\gamma\in\Gamma_{\leq}(\mu,\nu)$. Pick a Borel set
$A\subset\hat{\Omega}$, and suppose $\hat{\infty}\in A$. By definition
$\hat{\gamma}$ is a measure define on $\hat{\Omega}^{2}$, then we have

$\displaystyle\hat{\gamma}(A\times\hat{\Omega})$
$\displaystyle=\hat{\gamma}(A\setminus\\{\hat{\infty}\\}\times\Omega)+\hat{\gamma}(A\setminus\\{\hat{\infty}\\}\times\\{\hat{\infty}\\})+\hat{\gamma}(\\{\hat{\infty}\\}\times\Omega)+\hat{\gamma}(\\{\hat{\infty},\hat{\infty}\\})$
$\displaystyle=\gamma(A\setminus\\{\hat{\infty}\\}\times\Omega)+(\mu-(\pi_{1})_{\\#}\gamma)(A\setminus\\{\hat{\infty}\\})+(\nu-(\pi_{2})_{\\#}\gamma)(\Omega)+\gamma(\Omega^{2})$
$\displaystyle=(\pi_{1})_{\\#}\gamma(A\setminus\\{\hat{\infty}\\})+(\mu-(\pi_{1})_{\\#}\gamma)(A\setminus\\{\hat{\infty}\\})+\nu(\Omega)-(\pi_{2})_{\\#}\gamma(\Omega)+\gamma(\Omega^{2})$
$\displaystyle=\mu(A\setminus\\{\hat{\infty}\\})+\nu(\Omega)$
$\displaystyle=\hat{\mu}(A)$

Similarly, if $\hat{\infty}\notin A$, we have
$\hat{\gamma}(A\times\hat{\Omega})=\mu(A)=\hat{\mu}(A)$. Thus
$(\pi_{1})_{\\#}\hat{\gamma}=\hat{\mu}$ and similarly
$(\pi_{2})_{\\#}\hat{\gamma}=\hat{\nu}$. Therefore
$\hat{\gamma}\in\Gamma(\hat{\mu},\hat{\nu})$.

It is obvious that the mapping $T$ is injective since if
$\gamma_{1}\neq\gamma_{2}$ where
$\gamma_{1},\gamma_{2}\in\Gamma_{\leq}(\mu,\nu)$, then there exists one set
$B\subset\Omega^{2}$ such that $\gamma_{1}(B)\neq\gamma_{2}(B)$. Then
$\hat{\gamma}_{1}(B)=\gamma_{1}(B)\neq\gamma_{2}(B)=\hat{\gamma}_{2}(B)$.
Therefore, $\hat{\gamma}_{1}\neq\hat{\gamma}_{2}$.

Next we will show the surjectivity of $T$. Pick any
$\hat{\gamma}\in\Gamma(\hat{\mu},\hat{\nu})$, define $\gamma$ such that for
any $B\subset\Omega^{2}$, $\gamma(B)=\hat{\gamma}(B)$. We have
$\gamma\in\Gamma_{\leq}(\mu,\nu)$. Indeed, pick Borel set $A\subset\Omega$, we
have

$\gamma(A\times\Omega)=\hat{\gamma}(A\times\Omega)\leq\hat{\gamma}(A\times\hat{\Omega})=\hat{\mu}(A)=\mu(A).$

Thus $(\pi_{1})_{\\#}\gamma\leq\mu$, similarly we have
$(\pi_{2})_{\\#}\gamma\leq\nu$. Let $\hat{\gamma}_{1}=T(\gamma)$. We claim
$\hat{\gamma}=\hat{\gamma}_{1}$. Note, since
$\Omega^{2},\Omega\times\\{\hat{\infty}\\},\\{\hat{\infty}\\}\times\Omega,\\{\hat{\infty},\hat{\infty}\\}$
is a disjoint decomposition of $\hat{\Omega}^{2}$ (and all of them are
measurable), it is sufficient to prove $\hat{\gamma}(B)=\hat{\gamma}_{1}(B)$
for any Borel set $B$ which is a subset of these four sets.

Case 1: If $B\subset\Omega^{2}$, we have
$\hat{\gamma}_{1}(B)=\gamma(B)=\hat{\gamma}(B)$.

Case 2: If $B=A\times\\{\hat{\infty}\\}$ where $A\subset\Omega$ is Borel set,
then

$\displaystyle\hat{\gamma}(B)$
$\displaystyle=\hat{\gamma}(A\times\hat{\Omega})-\hat{\gamma}(A\times\Omega)$
$\displaystyle=\hat{\mu}(A)-\gamma(A\times\Omega)$
$\displaystyle=\mu(A)-(\pi_{1})_{\\#}\gamma(A)$
$\displaystyle=\hat{\gamma}_{1}(A\times\\{\hat{\infty}\\})=\hat{\gamma}_{1}(B)$

Similarly, if $B=\\{\hat{\infty}\\}\times A$ for some $A\subset\Omega$, we
have $\hat{\gamma}(B)=\hat{\gamma}_{1}(B)$.

Case 3: If $B=\\{(\hat{\infty},\hat{\infty})\\}$. Note, since
$\hat{\gamma}_{1}\in\Gamma(\tilde{\mu},\tilde{\nu})$ as we discussed above,
then $\hat{\gamma}(\hat{\Omega}^{2})=\hat{\gamma}_{1}(\hat{\Omega}^{2})$.
Additionally, by Cases 1 and 2 we have
$\hat{\gamma}(\Omega\times\Omega)=\hat{\gamma}_{1}(\Omega\times\Omega)$,
$\hat{\gamma}(\Omega\times\\{\hat{\infty}\\})=\hat{\gamma}_{1}(\Omega\times\\{\hat{\infty}\\})$
and
$\hat{\gamma}(\\{\hat{\infty}\\}\times\Omega)=\hat{\gamma}_{1}(\\{\hat{\infty}\\}\times\Omega)$.
Thus

$\displaystyle\hat{\gamma}(B)$
$\displaystyle=\hat{\gamma}(\hat{\Omega}^{2})-\hat{\gamma}(\Omega^{2})-\hat{\gamma}(\Omega\times\\{\hat{\infty}\\})-\hat{\gamma}(\\{\hat{\infty}\\}\times\Omega)$
$\displaystyle=\hat{\gamma}_{1}(\hat{\Omega}^{2})-\hat{\gamma}_{1}(\Omega^{2})-\hat{\gamma}_{1}(\Omega\times\\{\hat{\infty}\\})-\hat{\gamma}_{1}(\\{\hat{\infty}\\}\times\Omega)$
$\displaystyle=\hat{\gamma}_{1}(B)$ (13)

Hence, $\hat{\gamma}=\hat{\gamma}_{1}$ and thus that the mapping is
surjective.

Finally, we will show $\gamma$ is optimal in $\mathrm{OPT}_{\lambda}(\mu,\nu)$
if and only if $\hat{\gamma}$ is optimal in $\text{OT}(\hat{\mu},\hat{\nu})$
(defined in (11)). We let $C(\gamma),\hat{C}(\hat{\gamma})$ denote the
corresponding transportation cost of $\gamma,\hat{\gamma}$ with respect to the
OPT, OT problems, i.e.

$\hat{C}(\hat{\gamma})=\int\hat{c}(x,y)\,\mathrm{d}\gamma(x,y),\qquad
C(\gamma)=\int
c(x,y)\,\mathrm{d}\gamma(x,y)+\frac{\lambda}{2}\left(\mu(\Omega)-\pi_{1\\#}\gamma(\Omega)+\nu(\Omega)-\pi_{2\\#}\gamma(\Omega)\right).$
(14)

We have
$C(\gamma)=C(\hat{\gamma})+\frac{1}{2}\lambda(\mu(\Omega)+\nu(\Omega))$.
Combining with the fact the mapping is bijection, we have $\gamma$ is optimal
iff $\hat{\gamma}$ is optimal. ∎

##  Appendix C Correctness and complexity of Algorithms 1 and 2

###  C.1 Correctness

In this section we prove the correctness of Algorithms 1 and 2 as stated above
and we discuss how to deal with duplicate points. Extended versions of the
Algorithms with more sophisticated data structures and proper handling of
boundaries and duplicates are then given in Section C.2 together with a bound
on their worst-case complexity.

#### Preliminaries, induction strategy, cases 1 and 2.

Throughout this proof we are simply going to write $c_{i,j}$ for
$c(x_{i},y_{j})$. We assume that the point lists $\\{x_{i}\\}_{i=1}^{n}$ and
$\\{y_{j}\\}_{j=1}^{m}$ are sorted, but we now allow for duplicate points and
their handling will be addressed throughout this proof. Since $c(x,y)=h(x-y)$
for $h$ strictly convex, it is easy to verify that

$\displaystyle c_{i,j}+c_{k,l}\leq c_{i,l}+c_{k,j}$ (15)

if $i\leq k$ and $j\leq l$, with a strict inequality if $x_{i}<x_{k}$ and
$y_{j}<y_{l}$. This is known as Monge property [6]. The proof works via
induction in the iterations of the main loop of Algorithm 1. We will show that
prior to the iteration for $x_{k}$ / after completing the iteration for
$x_{k-1}$, the following holds:

  1. I.

$\Psi_{j}\leq\lambda$ for all $j\in[1:m]$.

  2. II.

For all $j\in[1:m]$, if $\Psi_{j}<\lambda$, then $y_{j}$ is assigned.

  3. III.

$\Phi_{i}\leq\lambda$ for all $i\in[1:n]$.

  4. IV.

For all $i\in[1:k-1]$, if $\Phi_{i}<\lambda$, then $x_{i}$ is currently
assigned.

  5. V.

All dual constraints $\Phi_{i}+\Psi_{j}\leq c_{i,j}$ for all $i\in[1:n]$,
$j\in[1:m]$ hold.

  6. VI.

For all $i\in[1:n]$, $j\in[1:m]$, whenever $x_{i}$ is assigned to $y_{j}$, one
has $\Phi_{i}+\Psi_{j}=c_{i,j}$.

  7. VII.

The assignment $L$ will be monotonous. I.e. if $L[i]\neq-1$,
$L[i^{\prime}]\neq-1$ for $i<i^{\prime}$, then $L[i]<L[i^{\prime}]$.

We initialize with $\Psi_{j}=\lambda$ for all $j$, $\Phi_{i}=-\infty$ for all
$i$, and empty assignment $L_{i}=-1$ for all $i$. Therefore, prior to the
first iteration, when $k=0$, all conditions are satisfied. Next, note that
items (I) and (II) will be satisfied throughout the algorithm, since entries
of $\Psi$ are only ever decreased during the algorithm; entries are only
decreased when the corresponding $y_{j}$ are assigned; and once a point
$y_{j}$ is assigned, it may become re-assigned, but it never becomes un-
assigned again. The claim that $y_{j}$ is never un-assigned is clear in all
cases apart from Case 3.1. In Case 3.1 it follows from property (VIII) below,
which implies that when $x_{i_{\Delta}}$ is un-assigned from, say,
$y_{j_{\Delta}}$, then $y_{j_{\Delta}}$ is re-assigned to $i_{\Delta}+1$
(since the assignment between $x_{i^{\prime}}$ and $y_{j^{\prime}}$ satisfies
$L[i^{\prime}]=j_{\min}+i^{\prime}-i_{\min}$, i.e. the assignment is
consecutive).

Throughout the algorithm let $j_{\textnormal{last}}$ be the largest index
among any assigned points $y_{j}$. We initially set $j_{\textnormal{last}}=-1$
when no $y_{j}$ is assigned. Since assigned $y_{j}$ do not get un-assigned
(merely re-assigned), $j_{\textnormal{last}}$ is non-decreasing.

######  Lemma C.1.

If $j_{\textnormal{last}}\neq-1$, then during the iteration of the main loop
of Algorithm 1, for any minimizer $j^{*}$ in line 3 one has $y_{j^{*}}\geq
y_{j_{\textnormal{last}}}$. In particular, $j^{*}$ can always be chosen such
that $j^{*}\geq j_{\textnormal{last}}$.

###### Proof.

If $j_{\textnormal{last}}=-1$ there is nothing to prove, since $j^{*}\geq 1$.
If $j_{\textnormal{last}}\neq-1$, then there must be some $i\in[1:k-1]$ such
that $L[i]=j_{\textnormal{last}}$ and therefore
$\Phi_{i}+\Psi_{j_{\textnormal{last}}}=c_{i,j_{\textnormal{last}}}$. After
adjusting $\Phi_{k}$ in line 4 one has $\Phi_{k}+\Psi_{j^{*}}=c_{k,j^{*}}$. By
dual feasibility we have in addition

$\displaystyle\Phi_{k}+\Psi_{j_{\textnormal{last}}}$ $\displaystyle\leq
c_{k,j_{\textnormal{last}}},$ $\displaystyle\Phi_{i}+\Psi_{j^{*}}$
$\displaystyle\leq c_{i,j^{*}}.$

Combining these four (in-)equalities we get

$\displaystyle c_{k,j^{*}}+c_{i,j_{\textnormal{last}}}\leq
c_{k,j_{\textnormal{last}}}+c_{i,j^{*}}.$

If $x_{k}>x_{i}$, then by (15) we have $y_{j^{*}}\geq
y_{j_{\textnormal{last}}}$. So $j^{*}<j_{\textnormal{last}}$ can only happen
if $y_{j^{*}}=y_{j_{\textnormal{last}}}$ and thus we may also choose
$j_{\textnormal{last}}$ as minimizing index. Therefore, we may impose the
constraint $j^{*}\geq j_{\textnormal{last}}$ in line 3. In the case
$x_{k}=x_{i}$, assume $j_{\textnormal{last}}$ would not be a minimal index in
line 3, i.e.

$c_{k,j^{*}}-\Psi_{j^{*}}<c_{k,j_{\textnormal{last}}}-\Psi_{j_{\textnormal{last}}}=\Phi_{i},$

where in the last equality we used $x_{k}=x_{i}$. Since
$L[i]=j_{\textnormal{last}}$ one must have that the dual constraint for
$(i,j_{\textnormal{last}})$ must be active. This would imply that the dual
constraint for $(i,j^{*})$ is violated, which contradicts the induction
hypothesis. ∎

Now during iteration $k$, the change of $\Phi_{k}$ in line 4, by construction,
preserves (III) and (V). Assume we enter Case 1. The assignment function $L$
is not changed, hence (VI) and (VII) remain preserved, and since
$\Phi_{k}=\lambda$, (IV) is extended to $i=k$. Assume we enter Case 2. Then we
have $L_{k}=j^{*}$, $\Phi_{k}+\Psi_{j^{*}}=c_{k,j^{*}}$ and
$\Phi_{k}<\lambda$. Hence, (VI) remains true, and (IV) is extended to $i=k$.
If we choose $j^{*}>j_{\textnormal{last}}$ (which is possible by Lemma C.1),
we preserve (VII).

#### Case 3.

We now turn to Case 3.

######  Lemma C.2.

In each iteration of the main loop of Algorithm 1, when we enter Case 3, i.e.
$\Phi_{k}<\lambda$ and $j^{*}=j_{\textnormal{last}}$, let $i$ be the index
such that $x_{i}$ is currently assigned to $y_{j_{\textnormal{last}}}$. Then
$x_{i}=x_{i^{\prime}}$ for all $i^{\prime}\in[i:k-1]$. If $i<k-1$, then one
must have $\Phi_{i}=\lambda$.

###### Proof.

Clearly, $i<k$ (since it was assigned during a previous iteration). In the
following, let
$f(x)=c(x,y_{j_{\textnormal{last}}})-\Psi_{j_{\textnormal{last}}}$, which is
convex in $x\in\mathbb{R}$. By (VI) we have
$\Phi_{i}=c(x_{i},y_{j_{\textnormal{last}}})-\Psi_{j_{\textnormal{last}}}=f(x_{i})\leq\lambda$,
by the current iteration of the main loop we have
$\Phi_{k}+\Psi_{j_{\textnormal{last}}}=c(x_{k},y_{j_{\textnormal{last}}})=f(x_{k})<\lambda$.
Let now $i^{\prime}\in[i+1,k-1]$. By monotonicity of $L$, (VII), if
$L[i^{\prime}]\neq-1$, then we would need
$L[i^{\prime}]>j_{\textnormal{last}}$, which contradicts the definition of
$j_{\textnormal{last}}$. Therefore $L[i^{\prime}]=-1$ and therefore by (IV) we
must have $\Phi_{i^{\prime}}=\lambda$. By (V) we must also have
$\Phi_{i^{\prime}}\leq f(x_{i^{\prime}})$, and by convexity of $f$,
$f(x_{i^{\prime}})$, since $f(x_{i})\leq\lambda$, $f(x_{k})<\lambda$, this can
only happen if $x_{i^{\prime}}=x_{i}$, and $\Phi_{i}=f(x_{i})=\lambda$. ∎

This means that if all points are distinct, then we must have $i=k-1$, and
find ourselves in the main loop of Algorithm 2, see below.

######  Remark C.3.

If points are not necessarily distinct (at least up to numerical rounding
errors), and we find $i<k-1$, then the situation can be remedied by setting
$L[i]=-1$, $L[k]=j^{*}$, which preserves (VI), (IV) and (VII). We will add
this to the algorithms in Section C.2.

We now study Algorithm 2 to resolve the conflict. In addition to the items
above, at the beginning of each iteration of the main loop of Algorithm 2 the
following is preserved:

  1. VIII.

There are indices $j_{\min}\leq j^{*}$, $i_{\min}\leq k-1$ with
$j^{*}-j_{\min}=(k-1)-i_{\min}$ such that $L[i_{\min}+r]=j_{\min}+r$ for
$r\in[0:(k-1)-i_{\min}]$,
$\Phi_{i+r}+\Psi_{j_{\min}+r-1}=c_{i+r,j_{\min}+r-1}$ for
$r\in[1:(k-1)-i_{\min}]$.

This is clearly true before the first iteration, when $i_{\min}=k-1$ and
$j_{\min}=j^{*}$. In each iteration of the main loop we then seek the largest
possible value $\Delta\geq 0$ such that by setting
$\Phi_{i}\leftarrow\Phi_{i}+\Delta$ for $i\in[i_{\min},k]$ and
$\Psi_{j}\leftarrow\Psi_{j}-\Delta$ for $j\in[j_{\min},j^{*}]$ we preserve all
items (I) to (VII). Clearly the delicate ones are (III) and (V). To preserve
the former, we ensure that $\Delta\leq\lambda_{\Delta}$. To preserve the
latter, we do not need to worry about the $\Psi_{j}$, $j\in[j_{\min}:j^{*}]$,
since they are decreased, but we need to consider all constraints
$\Phi_{i}+\Psi_{j}\leq c_{i,j}$ for $i\in[i_{\min}:k]$,
$j\in[1:j_{\min}-1]\cup[j^{*}+1:m]$. By the following lemma, we see that this
can be reduced to checking the two constraints for $(i_{\min},j_{\min}-1)$ and
$(k,j^{*}+1)$, which is the role of the variables $\alpha$ and $\beta$ in
Algorithm 2.

######  Lemma C.4.

In the above situation, one has that

$\displaystyle\min_{\begin{subarray}{c}i\in[i_{\min}:k],\\\
j\in[j^{*}+1:m]\end{subarray}}c_{i,j}-\Phi_{i}-\Psi_{j}$
$\displaystyle=c_{k,j^{*}+1}-\Phi_{k}-\Psi_{j^{*}+1},$
$\displaystyle\min_{\begin{subarray}{c}i\in[i_{\min}:k],\\\
j\in[1:j_{\min}-1]\end{subarray}}c_{i,j}-\Phi_{i}-\Psi_{j}$
$\displaystyle=c_{i_{\min},j_{\min}-1}-\Phi_{i_{\min}}-\Psi_{j_{\min}-1},$

if $j^{*}<m$ and $j_{\min}>1$ respectively.

###### Proof.

We start with the first equation and begin by showing that

$c_{i,j}-\Phi_{i}-\Psi_{j}\geq c_{k,j}-\Phi_{k}-\Psi_{j}$ (16)

for $i\in[i_{\min}:k]$, $j\in[j^{*}+1:m]$. We get this by combining
$\Phi_{i}\leq c_{i,j^{*}}-\Psi_{j^{*}}$, $\Psi_{j^{*}}=c_{k,j^{*}}-\Phi_{k}$
and the Monge property of cost matrix, (15), $c_{k,j}\leq
c_{i,j}+c_{k,j^{*}}-c_{i,j^{*}}$. Next, observe that $\Psi_{j}=\lambda$ for
$j\in[j^{*}+1:m]$, since these values have not yet been changed since the
initialization, and $\Psi_{j^{*}}\leq\lambda$. Also we know that
$\Phi_{k}=c_{k,j^{*}}-\psi_{j^{*}}\leq c_{k,j}-\psi_{j}$ for all $j\in[1:m]$.
Combining this, we get $c_{k,j^{*}}\leq c_{k,j}$ for $j\in[j^{*}+1:m]$. Since
$f:y\mapsto c(x_{k},y)$ is convex, and $f(y_{j^{*}})\leq f(y_{j})$ for $j\geq
j^{*}$, we must have that $f$ is non-decreasing after $y_{j^{*}}$, and
therefore among all indices $j\geq j^{*}$, the smallest one attains the
minimum.

Now we turn to the second equation. In complete analogy to (16) we show that
$c_{i,j}-\Phi_{i}-\Psi_{j}\geq c_{i_{\min},j}-\Phi_{i_{\min}}-\Psi_{j}$ for
$i\in[i_{\min}:k]$, $j\in[1:j_{\min}-1]$. Arguing as in Lemma C.1, we can show
a minimizing index $j$ can be chosen such that it is not smaller than
$\hat{j}$, where $\hat{j}$ is the largest index among the assigned $y_{j}$,
that is less than $j_{\min}$ (if such an assigned point exists, otherwise just
let $\hat{j}=0$ in the following). Consequently, all $y_{j}$ points in
$[\hat{j}+1:j_{\min}-1]$ must be unassigned and therefore have
$\Psi_{j}=\lambda$. Arguing then via the convexity of $c$ as in the previous
paragraph, we can show that a minimizing $j$ must be given by $j_{\min}-1$. ∎

The selection of Cases 3.1, 3.2 or 3.3 depends now on which of the three
bounds $\lambda_{\Delta}$, $\alpha$, or $\beta$ is smallest. Consequently,
each of the implied updates of the dual variables in the three cases preserves
the dual constraints and it is easy to see that by property (VIII) each of the
conflict resolutions in Cases 3.1, 3.2 and 3.3a preserve all other conditions
(I) to (VII). For instance, when $\lambda_{\Delta}$ is minimal, element
$x_{i_{\Delta}}$ becomes unassigned, however we then have
$\Phi_{i_{\Delta}}=\lambda$ as required by (IV).

We are left with discussing Case 3.3b, i.e. when $\beta$ is minimal among the
three bounds and $y_{j_{\min}-1}$ is already assigned. In complete analogy to
Lemma C.2 we can prove the following.

######  Lemma C.5.

In each iteration of the main loop of Algorithm 2, when we enter Case 3.3b,
i.e. $\Phi_{i_{\min}}<\lambda$ and $y_{j_{\min}-1}$ is already assigned, let
$i$ be the index such that $x_{i}$ is currently assigned to $y_{j_{\min}-1}$.
Then $x_{i}=x_{i^{\prime}}$ for all $i^{\prime}\in[i:i_{\min}-1]$. If
$i<i_{\min}-1$, then one must have $\Phi_{i}=\lambda$.

As above, this means that if all points are distinct, then $i=i_{\min}-1$, we
can set $i_{\min}\leftarrow i_{\min}-1$, $j_{\min}\leftarrow j_{\min}-1$, note
that we satisfy $\Phi_{i_{\min}}+\Psi_{j_{\min}}=c_{i_{\min},j_{\min}}$ and
thus preserve (VIII) before the next iteration in Algorithm 2.

######  Remark C.6.

If points are not all distinct and if $i<i_{\min}-1$, then we must have
$\Phi_{i}=\lambda$, thus we can unassign $x_{i}$ and $y_{j_{\min}-1}$, and
then proceed as if $y_{j_{\min}-1}$ were unassigned and resolve the conflict
as in Case 3.3a.

###  C.2 Full algorithm versions and complexity

We now give more complete pseudo code versions of the Algorithms 1 and 2, see
Algorithms 4 and 5. The main purpose is to reach a quadratic worst case time
complexity. Our algorithm can be seen as a specialization of the Hungarian
method, exploiting the particular one-dimensional structure of the cost and
dealing consistently with the option to discard mass for a cost $\lambda$. The
changes are explained below, subsequently some additional changes (for
duplicate and boundary handling) are described in plain text, and finally we
show how to determine the time complexity bound.

Input: $\\{x_{i}\\}_{i=1}^{n},\\{y_{j}\\}_{j=1}^{m},\lambda$

Output: $L$, $\Psi$, $\Phi$

1 Initialize $\Phi_{i}\leftarrow-\infty$ for $i\in[1:n]$,
$\Psi_{j}\leftarrow\lambda$ for $j\in[1:m]$ and $L_{i}\leftarrow-1$ for
$i\in[1:n]$,

2 $j_{\textnormal{last}}\leftarrow 1$

3 for  _$k=1,2,\ldots n$_ do

4
$j^{*}\leftarrow\operatorname{argmin}_{j\in[j_{\textnormal{last}}:m]}c(x_{k},y_{j})-\Psi_{j}$

5 $\Phi_{k}\leftarrow\min\\{c(x_{k},y_{j^{*}})-\Psi_{j^{*}},\lambda\\}$

6 if  _$\Phi_{k}=\lambda$_ then

7 [Case 1] No update on $L$

8 else if  _$j_{\min}-1$ unassigned_ then

9 [Case 2] $L_{k}\leftarrow j^{*}$, $j_{\textnormal{last}}\leftarrow j^{*}$

10 else

11 [Case 3] Run Algorithm 5.

12

Algorithm 4 opt-1d

Input: ($\\{x_{i}\\}_{i=1}^{n},\\{y_{j}\\}_{j=1}^{m}$, $k$, $j^{*}$,
$j_{\textnormal{last}}$, $L$, $\Phi$, $\Psi$)

Output: (Updated $L$, $\Phi$, $\Psi$, optimal for
$\textnormal{OPT}(\\{x_{i}\\}_{i=1}^{k},\\{y_{j}\\}_{j=1}^{m})$, and
$j_{\textnormal{last}}$)

1Initialize $i_{\min}\leftarrow k-1$, $j_{\min}\leftarrow j^{*}$.

2 Initialize $v\leftarrow 0$, $d_{j}\leftarrow 0$ for $j\in[1:m]$,
$d_{k}\leftarrow 0$, $d_{k-1}\leftarrow 0$.

3 $i_{\Delta}\leftarrow\operatorname{argmin}_{i\in[k-1:k]}(\lambda-\Phi_{i})$,
$\lambda_{\Delta}\leftarrow\lambda-\Phi_{i_{\Delta}}$

4 while  _true_ do

5 $\alpha\leftarrow c(x_{k},y_{j^{*}+1})-\Phi_{k}-v-\Psi_{j^{*}+1}$

6 $\beta\leftarrow
c(x_{i_{\min}},y_{j_{\min}-1})-\Phi_{i_{\min}}-\Psi_{j_{\min}-1}$

7 if  _$\lambda_{\Delta}\leq\min\\{\alpha,\beta\\}$_ then

8 [Case 3.1]

9 $v\leftarrow v+\lambda_{\Delta}$

10 for  _$i\in[i_{\min},k-1]$_ do

11 $\Phi_{i}\leftarrow\Phi_{i}+v-d_{i}$,
$\Psi_{L_{i}}\leftarrow\Psi_{L_{i}}-v+d_{i}$

12 $\Phi_{k}\leftarrow\Phi_{k}+v$

13 $L_{i_{\Delta}}\leftarrow-1$, $L_{k}\leftarrow j^{*}$

14 for  _$i\in[i_{\Delta}+1:k-1]$_ do

15 $L_{i}\leftarrow L_{i}-1$

16 return

17

18 else if  _$\alpha\leq\min\\{\lambda_{\Delta},\beta\\}$_ then

19 [Case 3.2]

20 $v\leftarrow v+\alpha$

21 for  _$i\in[i_{\min},k-1]$_ do

22 $\Phi_{i}\leftarrow\Phi_{i}+v-d_{i}$,
$\Psi_{L_{i}}\leftarrow\Psi_{L_{i}}-v+d_{i}$

23 $\Phi_{k}\leftarrow\Phi_{k}+v$

24 $L_{k}\leftarrow j^{*}+1$, $j_{\textnormal{last}}\leftarrow j^{*}+1$

25 return

26 else

27 $v\leftarrow v+\beta$

28 if  _$j_{\min}-1$ unassigned_ then

29 [Case 3.3a]

30 for  _$i\in[i_{\min},k-1]$_ do

31 $\Phi_{i}\leftarrow\Phi_{i}+v-d_{i}$,
$\Psi_{L_{i}}\leftarrow\Psi_{L_{i}}-v+d_{i}$

32 $\Phi_{k}\leftarrow\Phi_{k}+v$

33 $L_{i_{\min}}\leftarrow j_{\min}-1$, $L_{k}\leftarrow j^{*}$

34 for  _$i\in[i_{\min}+1:k-1]$_ do

35 $L_{i}\leftarrow L_{i}-1$

36 return

37

38 else

39 [Case 3.3b]

40 $d_{i_{\min}-1}\leftarrow v$,
$\lambda_{\Delta}\leftarrow\lambda_{\Delta}-\beta$,

41 $i_{\min}\leftarrow i_{\min}-1$, $j_{\min}\leftarrow j_{\min}-1$

42 if  _$\lambda-\Phi_{i_{\min}} <\lambda_{\Delta}$_ then

43 $\lambda_{\Delta}\leftarrow\lambda-\Phi_{i_{\min}}$, $i_{\Delta}\leftarrow
i_{\min}$

44

45

46

Algorithm 5 sub-opt-full

#### Implemented modifications compared to Algorithms 1 and 2.

Compared to Algorithm 1, in Algorithm 4 we have added the variable
$j_{\textnormal{last}}$ for improved handling of duplicate points (or limited
numerical precision), see Remark C.3. Note that initializing
$j_{\textnormal{last}}\leftarrow 1$, even when no points are yet assigned,
yields the desired behaviour. Additional adaptations related to this are
discussed in the paragraph below.

Compared to Algorithm 2, there are several adaptations to Algorithm 5.

The dual variables $\Phi$ and $\Psi$ are not updated during every loop of the
algorithm but only once, when the conflict is resolved. This is handled via
the auxiliary variable $v$ and the auxiliary array $d$. The former stores the
total increment that will need to be applied to $\Phi_{k}$ at the end, in
addition $d_{i}$ stores the value of $v$ at the time when $i$ was added to the
‘chain’, therefore $v-d_{i}$ will be the necessary increment of $\Phi_{i}$ at
the end. This trick (which is also known for the Hungarian method) removes the
necessity to loop over the whole chain to update the dual variables during
each iteration of the main loop in Algorithm 5 and thus reduces the worst case
time complexity from cubic to quadratic.

Similarly, the index of the dual variable $\Phi_{i}$ that is currently closest
to $\lambda$ is not determined from scratch during each iteration. Instead,
when case 3.3b is entered, the old best value is first reduced by $\beta$,
then compared with the new competitor $i_{\min}$ (after updating $i_{\min}$),
and updated if necessary.

#### Additional recommended modifications to algorithm.

In lines 5 and 5 of Algorithm 5 boundary checks should be added. E.g. $\alpha$
can only be set as described if $j^{*}<m$, otherwise it should be set to
$+\infty$. Likewise, $\beta$ can only be set as described if $j_{\min}>1$ and
should otherwise be set to $+\infty$. To keep track of which $y_{j}$ are
assigned one can use a boolean array of size $m$, initialized with false, and
entries corresponding to assigned points are set to true. This can be used to
distinguish between cases 3.3a and 3.3b. Alternatively, an ‘inverse’ version
of $L$ can be maintained, where $L^{-1}[j]$ will store the index $i$ of point
$x_{i}$ to which point $y_{j}$ is assigned (and $-1$ otherwise). This has to
be updated consistently with $L$. The latter will be useful when dealing with
duplicate points according to Remarks C.3 and C.6 at the beginnings of case 3
and case 3.3b respectively.

#### Worst case time complexity.

In Algorithm 4, initialization of the arrays $\Phi$, $\Psi$ and $L$ requires
$\Theta(n+m)$ steps. The main loop runs exactly $n$ times. Determining $j^{*}$
requires $O(m)$ steps (using the particular structure of $c$ and
$\Psi_{j}=\lambda$ for $j>j^{*}$ one could reduce this further, see Lemma C.4,
but we leave such optimizations for future work). Cases 1 and 2 take
$\Theta(1)$ steps. Let us now consider case 3 and Algorithm 5. Initialization
takes $\Theta(m)$ for setting up $d$ (although we note that this
initialization could be skipped). Cases 3.1, 3.2 and 3.3a are each entered
only once, right before termination of the sub-routine, and they have a
complexity of $O(n)$ (iterating over the chain for a fixed number of times to
adjust $L$ and the duals). Case 3.3b, as well as maintaining the variables
$\alpha$, $\beta$ and $\lambda_{\Delta}$ have a complexity of $\Theta(1)$ per
iteration and there are $O(n)$ iterations. Hence, Algorithm 5 in its current
form has a complexity of $O(\max\\{m,n\\})$, and therefore finally Algorithm 4
has a complexity of $O(n\cdot\max\\{m,n\\})$. During the proof we have pointed
out some potential for optimizing the algorithm for the regime when $n\ll m$.