# Anchor Space Optimal Transport: Accelerating Batch Processing of Multiple OT
Problems

Jianming Huang Graduate School of Fundamental Science and Engineering, WASEDA
University, 3-4-1 Okubo, Shinjuku-ku, Tokyo 169-8555, Japan (e-mail:
<EMAIL_ADDRESS> Xun Su Graduate School of Fundamental Science and
Engineering, WASEDA University, 3-4-1 Okubo, Shinjuku-ku, Tokyo 169-8555,
Japan (e-mail<EMAIL_ADDRESS>Zhongxi Fang Graduate School of
Fundamental Science and Engineering, WASEDA University, 3-4-1 Okubo, Shinjuku-
ku, Tokyo 169-8555, Japan (e-mail<EMAIL_ADDRESS>Hiroyuki Kasai
Department of Computer Science and Communication Engineering, WASEDA
University, 3-4-1 Okubo, Shinjuku-ku, Tokyo 169-8555, Japan (e-mail:
<EMAIL_ADDRESS>
###### Abstract

The optimal transport (OT) theory provides an effective way to compare
probability distributions on a defined metric space, but it suffers from cubic
computational complexity. Although the Sinkhorn’s algorithm greatly reduces
the computational complexity of OT solutions, the solutions of multiple OT
problems are still time-consuming and memory-comsuming in practice. However,
many works on the computational acceleration of OT are usually based on the
premise of a single OT problem, ignoring the potential common characteristics
of the distributions in a mini-batch. Therefore, we propose a translated OT
problem designated as the anchor space optimal transport (ASOT) problem, which
is specially designed for batch processing of multiple OT problem solutions.
For the proposed ASOT problem, the distributions will be mapped into a shared
anchor point space, which learns the potential common characteristics and thus
help accelerate OT batch processing. Based on the proposed ASOT, the
Wasserstein distance error to the original OT problem is proven to be bounded
by ground cost errors. Building upon this, we propose three methods to learn
an anchor space minimizing the distance error, each of which has its
application background. Numerical experiments on real-world datasets show that
our proposed methods can greatly reduce computational time while maintaining
reasonable approximation performance.

##  1 Introduction

As an important tool for comparing probability distributions, the optimal
transport (OT) theory [1, 2] provides a robust and relevant distance for
measuring two probability distributions on an already defined metric space,
which is so-called the Wasserstein distance. However, due to the high
computational complexity of $O(n^{3}\log(n))$ when solving the OT problem, it
is difficult to apply it to practical machine learning tasks. In recent years,
thanks to the computational improvements in solving OT problems with entropic
regularizers, computing an approximate Wasserstein distance with a negligible
error is brought down to a quadratic time complexity. Therefore, a lot of
research works with applications of the OT theory showed up in many domains of
machine learning, such as [3, 4], which introduced a mini-batch OT loss to the
generative models. In the graph learning domain, a framework of the
Wasserstein kernel was proposed, which utilizes the Wasserstein distance to
define distances or embeddings between graph-structured data [5, 6]. In the
sequence learning domain, [7] modified the original OT problem to generate a
solution with the characteristic of sequential data by adding temporal
regularization terms, while [8] proposed a method for sequence matching with
three auto-weighted OT subproblems. In the meantime, in order to compute the
solution to the OT problem more efficiently, many efforts have also been made.
This includes related works such as the sliced Wasserstein distance [9, 10],
which approximates the solution by solving one-dimensional OT problems on
projections of data.

There exist many machine learning tasks that need to solve multiple OT
problems as described in Definition 1 of Section 2. For example, in a graph
classification based on the Wasserstein graph kernels [5, 11, 12], we should
compute a pairwise Wasserstein graph distance matrix in order to instantiate a
precomputed kernel for the support vector machine. The pairwise Wasserstein
graph distance matrix contains distances of every pair of graphs of the
dataset, each of which is computed by solving a single OT problem. There are
other similar examples such as the $k$-NN sequence classification task based
on the order-preserving Wasserstein distance [7], and the distribution
clustering tasks based on the Wasserstein distance [13, 14]. However, many
works on the computational acceleration of OT are usually based on the premise
of a single OT problem, without considering the batch processing acceleration
of multiple OT problems. In this case, the conventional one-by-one
computation, as shown at the left part of Figure 1, to compute a pairwise
Wasserstein distance matrix is very inefficient because of the operation to
instantiate the cost matrix (C in Eq. (1)) for each distribution pair.

Although the Sinkhorn’s solver [2] provides a GPU-parallelizable solution to
accelerate the OT batch processing, such operation is only applicable to
distributions with fixed sizes. Under a situation where distributions are of
different sizes, the OT batch processing becomes more difficult. Fortunately,
this difficulty could be resolved by an improved Sinkhorn’s solver with a
block-diagonal stacking strategy as described in Section 2.4. Nevertheless, it
still requires more extra operations for block-diagonal stacking operations,
and most importantly, it still cannot avoid the huge time and space cost
caused by the cost matrix instantiation.

With considering the difficulties mentioned above, we assume that there exists
some potential common characteristics among the distributions of the given
dataset. With this assumption, we believe that conducting full-distribution
transportations is unnecessary in batch processing, which will make the
solution time-consuming and memory-consuming. Therefore, to resolve these
difficulties, we propose a translated OT problem, where the distributions are
transported based on an already learned and shared anchor point space.
Thereafter, the original OT problem is converted to a lightweight one between
a unique set of anchor points over the mapped distributions. As shown at the
right part of Figure 1, because all the translated OT problems can share the
same anchor point set, we do not need to instantiate the cost/kernel matrices
for each of them, which helps to reduce much computational time and required
memory for batch-processing. Moreover, because the distributions of the
translated OT problems have a fixed size whether they are originally size-
variable or not, it can be directly applicable to the Sinkhorn solver with GPU
parallelization, without the need of the block-diagonal stacking. We name it
the anchor space optimal transport (ASOT) problem. Furthermore, we derived and
proved the upper bound of the Wasserstein distance error of our proposed ASOT
problem to the original OT problem. We also proposed a metric learning
framework (ASOT-ML) to learn the anchor space that minimizes the upper bound
of Wasserstein distance error. However, we found that ASOT-ML has several
limitations, which make it difficult to use under some conditions. Therefore,
we further propose a lightweight version based on the $k$-means clustering
(ASOT-$k$), and a deep dictionary learning version (ASOT-DL).

While there are similar works mapping distributions into anchor spaces or
subspaces [15, 16], these focus on enhancing OT’s robustness rather than
considering multiple OT problems. For example, the latent OT [15] divides the
OT problem into three subproblems of OT: source space to source anchor space,
source anchor space to target anchor space, and target anchor space to target
space, which costs much greater computational time than OT. Conversely, our
proposed ASOT problem only contains one OT problem between the learned anchor
points.

Our contributions are threefold: (i) We propose an anchor space optimal
transport (ASOT) problem, which is translated from an optimal transport (OT)
problem but benefits from no need of pairwise cost/kernel instantiation and
the always fixed distribution size. The proposed ASOT greatly reduces the
computational costs of the batch processing of multiple OT problems. (ii) We
theoretically derived and proved the upper bound of the Wasserstein distance
error between our ASOT and OT, which makes our approximation errors
predictable. Moreover, to minimize such an upper bound, we propose three
methods for anchor space learning (ASOT-ML, ASOT-$k$, and ASOT-DL) for dealing
with different situations. (iii) We evaluate distance approximation errors and
computational time by conducting numerical experiments on several widely-used
real-world datasets. These results show that our proposed approachs have a
great reduction of computational time with reasonable approximation
performances.

Figure 1: Illustration of the two strategies to solve multiple OT problems:
the one-by-one computation (left) and the batch processing with the proposed
ASOT problem (right). Assume a multiple OT problem set
$\mathcal{M}(\\{\mathcal{D}_{1},\mathcal{D}_{2},\mathcal{D}_{3}\\},\\{(1,2),(1,3),(2,3)\\})$,
let $\mbox{\bf T}_{i,j}$ denotes the transportation plan from
$\mathcal{D}_{i}$ to $\mathcal{D}_{j}$. To solves the multiple problems, the
conventional one-by-one computation needs the instantiation of cost matrices
($\mbox{\bf C}_{1,2},\mbox{\bf C}_{1,3},\mbox{\bf C}_{2,3}$) for each
distribution pairs. On the other hand, our proposed method learns an anchor
space for all the distributions. Thereafter, the original distributions are
mapped into the learned anchor space as
$\mathcal{D}_{1}^{\prime},\mathcal{D}_{2}^{\prime},\mathcal{D}_{3}^{\prime}$,
which is a probability simplex based on a learnable anchor point set
$\mathcal{W}$. When conducting batch processing with ASOT, all the translated
OT problems share the same cost matrix $\mbox{\bf C}_{s}$, which is a pairwise
distance matrix of the anchor points in $\mathcal{W}$. Therefore, there is no
need to instantiate cost matrices for each distribution pairs.

##  2 Preliminaries

We use lowercase letters to represent scalars, such as $a,b$, and $c$.
Lowercase letters with bold typeface such as
$\mbox{\boldmath$a$},\mbox{\boldmath$b$},\mbox{\boldmath$c$}$ represent the
vectors, uppercase letters with bold typeface such as $\mbox{\bf A},\mbox{\bf
B},\mbox{\bf C}$ represent the matrices. For the needs of slicing matrices, we
represent the element of the $i$-th row and $j$-th columns of A as $\mbox{\bf
A}(i,j)$. Furthermore, $\mbox{\bf A}(i,:)$ denotes the vector of $i$-th row of
A, similarly, $\mbox{\bf A}(:,j)$ denotes the vector of $j$-th column.
$\mbox{\boldmath$a$}(i)$ denotes the $i$-th element of the vector $a$.
Uppercase letters with calligraphy typeface denote the set, such as
$\mathcal{A},\mathcal{B},\mathcal{C}$, $|\mathcal{A}|$ denotes the size of
$\mathcal{A}$. $\mathbb{R}$ denotes the real number set, $\mathbb{R}_{+}$
denotes the non-negative real number set. $\mbox{\boldmath$1$}_{n}$ denotes an
all-one vector. $\langle\mbox{\bf A},\mbox{\bf B}\rangle=\sum_{i,j}\mbox{\bf
A}(i,j)\mbox{\bf B}(i,j)$ denotes the Frobenius dot product.
$\Delta^{p}_{n}:=\\{\mbox{\boldmath$a$}\in\mathbb{R}^{n}:a_{i}\geq
0,\sum_{i=1}^{n}a_{i}=1\\}$ denotes the $n$-dimensional probability simplex,
$\Delta^{u}_{n}:=\\{\mbox{\boldmath$a$}\in\mathbb{R}^{n}:a_{i}\geq
0,\sum_{i=1}^{n}a_{i}\leq 1\\}$ denotes the $n$-dimensional unit simplex.
$[N]$ denotes the integer set $\\{1,2,\cdots,N\\}$.

###  2.1 Optimal transport (OT) problem

The optimal transport (OT) problem [1, 2] aims to find an OT plan between two
probability measures defined on a given metric space. For the main topic
discussed in this paper, we only consider the OT problem based on discrete
probability distributions, which could be viewed as finite sets of sampling
points with feature vectors. Assume that two discrete probability
distributions on different simplexes are given as
$\mbox{\boldmath$a$}\in\Delta_{n}^{p}$ and
$\mbox{\boldmath$b$}\in\Delta_{m}^{p}$. Assume that the sampling points of
$\mbox{\boldmath$a$},\mbox{\boldmath$b$}$ are on the same metric space, where
the distance between the sampling points could be measured by a given distance
metric. Let $\mbox{\bf C}\in\mathbb{R}^{n\times m}$ be the ground cost matrix
where $\mbox{\bf C}_{i,j}$ denotes the distance measured between the $i$-th
sampling point of $a$ and $j$-th one of $b$. The OT problem between $a$ and
$b$ is defined as

$\mathop{\rm minimize}_{{\bf
T}\in\mathcal{U}(\bm{a},\bm{b})}~{}\langle\mbox{\bf T},\mbox{\bf C}\rangle,$
(1)

where $\mathcal{U}(\mbox{\boldmath$a$},\mbox{\boldmath$b$})=\\{\mbox{\bf
T}\in\mathbb{R}_{+}^{n\times m}:\mbox{\bf
T}\mbox{\boldmath$1$}_{m}=\mbox{\boldmath$a$},\mbox{\bf
T}^{T}\mbox{\boldmath$1$}_{n}=\mbox{\boldmath$b$}\\}$ denotes the set of all
transportation plans between $a$ and $b$. The OT problem derives a new
distance measuring the optimal transportation cost named Wasserstein distance.
The $p$-Wasserstein distance between two discrete probability distributions
$\mbox{\boldmath$a$},\mbox{\boldmath$b$}$ is defined as
$W_{p}(\mbox{\boldmath$a$},\mbox{\boldmath$b$},\mbox{\bf C}):=\mathop{\rm
minimize}_{{\bf T}\in\mathcal{U}(\bm{a},\bm{b})}~{}\langle\mbox{\bf
T},\mbox{\bf C}^{p}\rangle^{\frac{1}{p}}$. A wide application of the
Wasserstein distance in specific machine learning tasks is to measure the
1-Wasserstein distance between two feature (sample) sets. Given two sets of
$d$-dimensional features vectors
$\mathcal{X}=\\{\mbox{\boldmath$x$}_{i}\\}_{i=1}^{n},\mathcal{Y}=\\{\mbox{\boldmath$y$}_{j}\\}_{j=1}^{m}$,
where $\mbox{\boldmath$x$}_{i},\mbox{\boldmath$y$}_{j}\in\mathbb{R}^{d}$,
$|\mathcal{X}|=n,|\mathcal{Y}|=m$ and their probability (mass) vector
$\mbox{\boldmath$a$}\in\Delta_{n}^{p}$ and
$\mbox{\boldmath$b$}\in\Delta_{m}^{p}$, respectively. Let
$\mathcal{D}_{x}=(\mathcal{X},\mbox{\boldmath$a$}),\mathcal{D}_{y}=(\mathcal{Y},\mbox{\boldmath$b$})$
denotes the two distributions based on these two sample sets. To measure the
optimal transportation cost between $\mathcal{D}_{x}$ and $\mathcal{D}_{y}$,
then the 1-Wasserstein distance between $\mathcal{D}_{x}$ and
$\mathcal{D}_{y}$ is denoted as

$W_{1}(\mathcal{D}_{x},\mathcal{D}_{y},d_{\rm S}):=\mathop{\rm minimize}_{{\bf
T}\in\mathcal{U}\left(\bm{a},\bm{b}\right)}\langle\mbox{\bf T},\mbox{\bf
C}\rangle,$ (2)

where $d_{\rm S}:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}_{+}$ denotes
the given ground cost metric, $\mbox{\bf C}(i,j):=d_{\rm
S}(\mbox{\boldmath$x$}_{i},\mbox{\boldmath$y$}_{j})$. It should be noted that,
because in a practical situation, there is usually no prior information about
the probability distributions that are already known, there is a common way to
set the probability vectors as uniform distributions as
$\mbox{\boldmath$a$}=\frac{\bm{1}_{n}}{n},\mbox{\boldmath$b$}=\frac{\bm{1}_{m}}{m}$.

###  2.2 Sinkhorn’s algorithm

Because the OT problem expressed in Eq. (1) could be represented in the
canonical form of a linear programming problem, which is a convex optimization
problem, it can be solved in the complexity of $O(n^{3}\log(n))$ in the case
of $n=m$. Another more popular scheme for solving OT problems is Sinkhorn’s
algorithm. Adding an entropic regularization term $H(\mbox{\bf
T}):=-\sum_{i,j}\mbox{\bf T}(i,j)(\log(\mbox{\bf T}(i,j))-1)$ to Eq. (1)
obtains an entropically regularized OT (eOT) problem as

$W^{e}_{1}(\mathcal{D}_{x},\mathcal{D}_{y}):=\mathop{\rm minimize}_{{\bf
T}\in\mathcal{U}(\bm{a},\bm{b})}~{}\langle\mbox{\bf T},\mbox{\bf
C}\rangle-\varepsilon H(\mbox{\bf T}),$

where $\varepsilon>0$ is a regularization paramter. This can be solved by
introducing Sinkhorn’s fixed-point iterations in complexity of
$O(n^{2}\log(n)\tau^{-3})$ with setting $\varepsilon=\frac{4\log(n)}{\tau}$
[17, 18, 19, 20]. In Sinkhorn’s algorithm, the transportation plan could be
decomposed as $\mbox{\bf T}=\mbox{\boldmath$u$}^{T}\mbox{\bf
K}{\mbox{\boldmath$v$}}$, where $\mbox{\bf K}:=\exp(-\frac{{\bf
C}}{\varepsilon})\in\mathbb{R}^{n\times m}$ is a Gibbs kernel computed from C
and
$\mbox{\boldmath$u$}\in\mathbb{R}^{n}_{+},\mbox{\boldmath$v$}\in\mathbb{R}^{m}_{+}$
are two scaling variables to be updated. Then the updating equations of the
$j$-th iteration of the Sinkhorn’s algorithm are defined as
$\mbox{\boldmath$u$}^{(j+1)}=\frac{\mbox{\boldmath$a$}}{\mbox{\bf
K}\mbox{\boldmath$v$}^{(j)}}$ and
$\mbox{\boldmath$v$}^{(j+1)}=\frac{\mbox{\boldmath$b$}}{\mbox{\bf
K}^{T}\mbox{\boldmath$u$}^{(j+1)}}$.

###  2.3 Multiple optimal transport problem set

We present our definition of the multiple OT problem set, which is the main
application scenario of the proposal.

######  Definition 1 (Multiple optimal transport problem set).

Assume a discrete distribution set of size $N$ as
$\mathbb{D}:=\\{\mathcal{D}_{1},\cdots,\mathcal{D}_{N}\\}$. Given a pair set
$\mathcal{P}\subset\\{(i,j):i,j\in[N],i\neq j\\}$, the multiple OT problem set
is defined as

$\mathcal{M}(\mathbb{D},\mathcal{P}):=\left\\{W_{1}(\mathcal{D}_{i},\mathcal{D}_{j},d_{\rm
S}):(i,j)\in\mathcal{P};\mathcal{D}_{i},\mathcal{D}_{j}\in\mathbb{D}\right\\},$
(3)

where $d_{\rm S}:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}_{+}$ denotes
the ground cost metric function. $d$ denotes the number of dimensions of each
sample.

Problems of the multiple OT problem set can be viewed as single OT problems of
the given distribution set $\mathbb{D}$. For the cases where distributions in
$\mathbb{D}$ are of the same size, such problems can be solved parallelly with
the GPU-parallelizable Sinkhorn’s algorithm. When distribution sizes are
allowed to vary, we cannot directly utilize the GPU-parallelizable Sinkhorn’s
algorithm and should instead utilize the block-diagonal stacking strategy.

###  2.4 The block-diagonal stacking strategy

Under the condition of size-fixed distributions, the updating equations of the
Sinkhorn’s algorithm can be easily batched with stacking
$\mbox{\boldmath$u$},\mbox{\boldmath$v$},\mbox{\boldmath$a$},\mbox{\boldmath$b$}$
and K along a new batch dimension to create a volume. Although this strategy
is inapplicable to size-variable distributions, we can adopt a block-diagonal
stacking strategy where the kernel matrices K are stacked along the diagonal,
while other vectors are stacked vertically instead of creating a new
dimension. Thereafter, substituting the batched matrices and vectors into the
original equation obtains a matrix of batched solutions. Such batch processing
is also efficiently GPU-parallelizable with the release of dedicated CUDA
kernels [21].

##  3 Anchor Space Optimal Transport (ASOT)

This section introduces the definition of our proposed ASOT problem. We then
present the theorem of upper bound of absolute 1-Wasserstein distance error
between ASOT and normal OT.

###  3.1 Motivation and definitions

As described in the previous sections, there are two weak points of the
conventional one-by-one computation that we focus on: (i) the need of cost
matrices instantiation for every distribution pairs; (ii) the need of extra
operations to parallelize with size-variable distributions. To resolve these
two difficulties, we believe that we can take advantages of some potential
common characterustics of the given distributions in the multiple OT problem
set. In order to remove the cost matrices instantiation described in (i), we
design a so-called ASOT problem to approximate the original OT problem. It can
be viewed as a new one with a low-rank representation of the original cost
matrix. The most improtant point is that it has a fixed and shareable
component (Lemma 1). For the difficulty of (ii), in our proposed ASOT problem,
the transportation is only based on the distributions of the learned anchor
point set, which have the fixed size equal to the number of anchor points,
thus do not require the operation of block-diagonal stacking. In summary, we
propose to learn a so-called anchor space for a batch of distributions that
allow variable size, over which our defined ASOT problems could be solved
efficiently. Meanwhile, it also has an upper bounded cost error to the
original OT problem. We present the relative definitions of ASOT as follows:

######  Definition 2 (Anchor space).

Assume a metric space on $\mathbb{R}^{d}$ with a metric function $d_{\rm
AS}:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}_{+}$, denoted as
$S(\mathbb{R}^{d},d_{\rm AS})$. Let
$\mathcal{W}=\\{\mbox{\boldmath$w$}_{i}\\}_{i=1}^{k}$ be a set of $k$ points
on $S(\mathbb{R}^{d},d_{\rm AS})$, where
$\mbox{\boldmath$w$}_{i}\in\mathbb{R}^{d}$. An anchor space of $\mathcal{W}$
is defined as a $k$-dimensional probability simplex, denoted as $S_{\rm
A}(\mathcal{W})=\Delta_{k}^{p}$. The points of $\mathcal{W}$ are called the
anchor points.

######  Definition 3 (Anchor space optimal transport (ASOT) problem).

Let $d_{\rm S},d_{\rm AS}:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}_{+}$
be two metric functions on $\mathbb{R}^{d}$. Assume two discrete distributions
$\mathcal{D}_{x}=(\mathcal{X},\mbox{\boldmath$a$}),\mathcal{D}_{y}=(\mathcal{Y},\mbox{\boldmath$b$})$
on a source metric space $S(\mathbb{R}^{d},d_{\rm S})$, where
$\mathcal{X}=\\{\mbox{\boldmath$x$}_{i}\\}_{i=1}^{n},\mathcal{Y}=\\{\mbox{\boldmath$y$}_{j}\\}_{j=1}^{m}$
denotes the sample sets and
$\mbox{\boldmath$x$}_{i},\mbox{\boldmath$y$}_{j}\in\mathbb{R}^{d}$,
$\mbox{\boldmath$a$}\in\Delta^{p}_{n},\mbox{\boldmath$b$}\in\Delta^{p}_{m}$
denotes their probability (mass) vectors, respectively. Given a set of $k$
anchor points $\mathcal{W}=\\{\mbox{\boldmath$w$}_{i}\\}_{i=1}^{k}$ on a new
metric space $S(\mathbb{R}^{d},d_{\rm AS})$, where
$\mbox{\boldmath$w$}_{i}\in\mathbb{R}^{d}$. And given a mapping function
$P:\mathbb{R}^{d}\to S_{\rm A}(\mathcal{W})$. Then, we define the
$1$-Wasserstein distance of ASOT problem over the anchor space $S_{\rm
A}(\mathcal{W})$ as

$\displaystyle W_{\rm AS}(\mathcal{D}_{x},\mathcal{D}_{y},\mathcal{W},P,d_{\rm
AS})\ :=\ \mathop{\rm minimize}_{{\bf
T}_{s}\in\mathcal{U}(\bm{a}^{\prime},\bm{b}^{\prime})}~{}\langle\mbox{\bf
T}_{s},\mbox{\bf C}_{s}\rangle,$ (4)
$\displaystyle\mbox{\boldmath$a$}^{\prime}=\frac{1}{n}\sum_{i=1}^{n}P(\mbox{\boldmath$x$}_{i})\mbox{\boldmath$a$}(i),\quad\mbox{\boldmath$b$}^{\prime}=\frac{1}{m}\sum_{j=1}^{m}P(\mbox{\boldmath$y$}_{j})\mbox{\boldmath$b$}(j),$
(5)

where $W_{\rm AS}(\mathcal{D}_{x},\mathcal{D}_{y},\mathcal{W},P,d_{\rm AS})$
denotes the $1$-Wasserstein distance of the ASOT problem, and $\mbox{\bf
T}_{s}\in\mathbb{R}^{k\times k}$ denotes the transportation plan between the
anchor points of $\mathcal{W}$. $\mbox{\bf C}_{s}\in\mathbb{R}^{k\times k}$
denotes the pairwise cost matrix of the anchor points of $\mathcal{W}$, where
$\mbox{\bf C}_{s}(i,j):=d_{\rm
AS}(\mbox{\boldmath$w$}_{i},\mbox{\boldmath$w$}_{j})$.

We further define an entropic variant of the ASOT problem envisioning the
application of Sinkhorn’s algorithm for GPU parallelization.

######  Definition 4 (Entropic anchor space optimal transport (eASOT)
problem).

Follow the same definitions of Definition 3. Then, an entropic variant problem
of Eq. (4) is defined as

$W^{e}_{\rm AS}(\mathcal{D}_{x},\mathcal{D}_{y},\mathcal{W},P,d_{\rm
AS}):=\mathop{\rm minimize}_{{\bf
T}_{s}\in\mathcal{U}(\bm{a},\bm{b})}~{}\langle\mbox{\bf T}_{s},\mbox{\bf
C}_{s}\rangle-\varepsilon H(\mbox{\bf T}_{s}).$ (6)

This is called the entropic ASOT problem denoted as eASOT.

###  3.2 Properties of ASOT

In this subsection, we show several important and useful properties of our
proposed ASOT.

ASOT meets the definition of an OT problem. Reviewing the definition of an OT
problem, it is obvious that ASOT is also an OT problem over two new
distributions
$\mathcal{D}_{x}^{\prime}=(\mathcal{W},\mbox{\boldmath$a$}^{\prime}),\mathcal{D}_{y}^{\prime}=(\mathcal{W},\mbox{\boldmath$b$}^{\prime})$
with a ground cost matrix $\mbox{\bf C}_{s}$ computed by a distance measure
$d_{\rm AS}$. This property allows us to solve an ASOT problem in the same way
as OT does. The minimum cost of the ASOT problem is also guaranteed to have
the same properties as the $1$-Wasserstein distance.

ASOT always has fixed distribution sizes. As shown in the equations of
Definition 3, the sizes of distributions
$\mathcal{D}_{x}^{\prime},\mathcal{D}_{y}^{\prime}$, the transportation plan
$\mbox{\bf T}_{s}$ and the cost matrix $\mbox{\bf C}_{s}$ are always fixed,
which only depends on the anchor space $S_{\rm A}(\mathcal{W})$. Therefore,
this property makes the eASOT problem efficiently solvable with GPU
parallelization.

Upper bounded error to OT. In order to translate an arbitrary OT problem into
the form of our proposed ASOT problem with ensuring a controllable cost error.
We present the upper bound of their absolute 1-Wasserstein distance error as
the following Proposition 1. Before that, we first prepare the following Lemma
1, which also reveals another important property of our proposed ASOT problem:
its equivalence to a specific OT problem of uniform distributions.

######  Lemma 1 (Equivalence lemma).

The following OT problem is equivalent to the ASOT problem $W_{\rm
AS}(\mathcal{D}_{x},\mathcal{D}_{y},\mathcal{W},P,d_{\rm AS})$, where
$\widehat{\mbox{\bf C}}=\mbox{\bf Z}_{x}\mbox{\bf C}_{s}\mbox{\bf Z}_{y}^{T}$.

$\mathop{\rm minimize}_{{\bf T}\in\mathcal{U}(\bm{a},\bm{b})}\langle\mbox{\bf
T},\widehat{\mbox{\bf C}}\rangle,$ (7)

where
$\mathcal{D}_{x}=(\mathcal{X},\mbox{\boldmath$a$}),\mathcal{D}_{y}=(\mathcal{Y},\mbox{\boldmath$b$})$.

###### Proof.

Let $W_{T}:=\langle\mbox{\bf T},\widehat{\mbox{\bf C}}\rangle$ in Eq. (7). It
can be rewritten as following equation with substituting $\widehat{\mbox{\bf
C}}=\mbox{\bf Z}_{x}\mbox{\bf C}_{s}\mbox{\bf Z}_{y}^{T}$ into it.

$W_{T}=\sum_{u,v\in[k]}\sum_{i\in[n],j\in[m]}\mbox{\bf
Z}_{x}(i,u)\cdot\mbox{\bf Z}_{y}(j,v)\cdot\mbox{\bf T}(i,j)\cdot\mbox{\bf
C}_{s}(u,v).$

Let $\mbox{\bf T}_{s}=\mbox{\bf Z}^{T}_{x}\mbox{\bf T}\mbox{\bf
Z}_{y}\in\mathbb{R}^{k\times k}$, where $k:=|\mathcal{W}|$ denotes the number
of anchor points. Substituting it into the above equation obtains

$W_{T}=\sum_{u,v\in[k]}\sum_{i,j\in[k]}\mbox{\bf T}_{s}(i,j)\cdot\mbox{\bf
C}_{s}(u,v)=\langle\mbox{\bf T}_{s},\mbox{\bf C}_{s}\rangle.$

Substituting it into Eq. (7) obtains $\mathop{\rm minimize}_{{\bf
T}\in\mathcal{U}(\bm{a},\bm{b})}\langle\mbox{\bf T}_{s},\mbox{\bf
C}_{s}\rangle$. To reparameterize this problem with $\mbox{\bf T}_{s}$, we
first right multiply it with $\mbox{\boldmath$1$}_{k}$. Then we have
$\mbox{\bf T}_{s}\mbox{\boldmath$1$}_{k}=\mbox{\bf Z}^{T}_{x}\mbox{\bf
T}\mbox{\bf Z}_{y}\mbox{\boldmath$1$}_{k}$. Furthermore, because each row of
$\mbox{\bf Z}_{x}$ and $\mbox{\bf Z}_{y}$ is a simplex vector, we have
$\mbox{\bf Z}_{y}\mbox{\boldmath$1$}_{k}=\mbox{\boldmath$1$}_{m}$. Then we
have $\mbox{\bf T}_{s}\mbox{\boldmath$1$}_{k}=\mbox{\bf Z}_{x}^{T}\mbox{\bf
T}\mbox{\boldmath$1$}_{m}=\frac{1}{n}\mbox{\bf
Z}^{T}_{x}\mbox{\boldmath$a$}=\mbox{\boldmath$a$}^{\prime}$, where
$\mbox{\boldmath$a$}^{\prime}$ denotes the same one in Definition 3. In the
same way, we can obtain $\mbox{\bf
T}_{s}^{T}\mbox{\boldmath$1$}_{k}=\mbox{\boldmath$b$}^{\prime}$, where
$\mbox{\boldmath$b$}^{\prime}$ denotes the same one in Definition 3. Finally,
Eq. (7) can be rewritten as follows.

$\mathop{\rm minimize}_{{\bf
T}_{s}\in\mathcal{U}(\bm{a}^{\prime},\bm{b}^{\prime})}\langle\mbox{\bf
T}_{s},\mbox{\bf C}_{s}\rangle.$

∎

######  Proposition 1 (Upper bound of the absolute 1-Wasserstein distance
error).

Assume an OT problem for 1-Wasserstein distance
$W_{1}(\mathcal{D}_{x},\mathcal{D}_{y},d_{\rm S})$, with $\mbox{\bf
C}\in\mathbb{R}^{n\times m}$ as the ground cost matrix and $d_{\rm S}$ as the
metric function. Define an ASOT problem $W_{\rm
AS}(\mathcal{D}_{x},\mathcal{D}_{y},\mathcal{W},P,d_{\rm AS})$ of Definition 3
with given $\mathcal{W},P$ and $d_{\rm AS}$. Let $\mbox{\bf
Z}_{x}\in\mathbb{R}^{n\times k},\mbox{\bf Z}_{y}\in\mathbb{R}^{m\times k}$ be
the matrices where $\mbox{\bf Z}_{x}(i,:)=P(\mbox{\boldmath$x$}_{i}),\mbox{\bf
Z}_{y}(j,:)=P(\mbox{\boldmath$y$}_{j})$. We define a reconstructed cost matrix
$\widehat{\mbox{\bf C}}=\mbox{\bf Z}_{x}\mbox{\bf C}_{s}\mbox{\bf
Z}_{y}^{T}\in\mathbb{R}^{n\times m}$. Then the following inequality holds.

$|W_{\rm AS}(\mathcal{D}_{x},\mathcal{D}_{y},\mathcal{W},P,d_{\rm
AS})-W_{1}(\mathcal{D}_{x},\mathcal{D}_{y},d_{\rm S})|<\|\widehat{\mbox{\bf
C}}-\mbox{\bf C}\|_{1}.$ (8)

###### Proof.

Let $\mbox{\bf T}^{\star},\widehat{\mbox{\bf T}}^{\star}$ be the optimal T of
Eqs. (2) and (7), respectively. According to Lemma 1, we have $|W_{\rm
AS}-W_{1}|=|\sum_{i,j}\widehat{\mbox{\bf T}}^{\star}(i,j)\widehat{\mbox{\bf
C}}(i,j)-\sum_{i,j}\mbox{\bf T}^{\star}(i,j)\mbox{\bf C}(i,j)|$. In the case
of $W_{\rm AS}>W_{1}^{u}$, we have the following inequalities.

$\displaystyle|W_{\rm AS}-W_{1}|$ $\displaystyle=$
$\displaystyle\sum_{i,j}\widehat{\mbox{\bf T}}^{\star}(i,j)\widehat{\mbox{\bf
C}}(i,j)-\sum_{i,j}\mbox{\bf T}^{\star}(i,j)\mbox{\bf C}(i,j)$
$\displaystyle\leq$ $\displaystyle\sum_{i,j}\mbox{\bf
T}^{\star}(i,j)(\widehat{\mbox{\bf C}}(i,j)-\mbox{\bf C}(i,j))$
$\displaystyle\leq$ $\displaystyle\sum_{i,j}\mbox{\bf
T}^{\star}(i,j)|\widehat{\mbox{\bf C}}(i,j)-\mbox{\bf C}(i,j)|$
$\displaystyle<$ $\displaystyle\|\widehat{\mbox{\bf C}}-\mbox{\bf
C}\|_{\infty}<\|\widehat{\mbox{\bf C}}-\mbox{\bf C}\|_{1}.$

In the case of $W_{\rm AS}\leq W_{1}$, we have $|W_{\rm
AS}-W_{1}|=|\sum_{i,j}\mbox{\bf T}^{\star}(i,j)\mbox{\bf
C}(i,j)-\sum_{i,j}\widehat{\mbox{\bf T}}^{\star}(i,j)\widehat{\mbox{\bf
C}}(i,j)|$. Replacing $\mbox{\bf T}^{\star}(i,j)$ in the 2nd and the 3rd rows
above with $\widehat{\mbox{\bf T}}^{\star}(i,j)$ can obtain the desired
result. ∎

The above two proofs are both proof sketchs, while the complete proofs of
Lemma 1 and Proposition 1 are included in the supplementary file. Based on
Proposition 1, we can translate a given OT problem into an ASOT problem by
minimizing the upper bound $\|\widehat{\mbox{\bf C}}-\mbox{\bf C}\|_{1}$,
which can be viewed as a function of three variables $\mathcal{W},P$ and
$d_{\rm AS}$. To solve this optimization problem, we propose three methods for
different application backgrounds, one of which is based on the metric
learning and the other two are based on $k$-means clustering and the deep
dictionary learning, which will be discussed in Section 4.

###  3.3 GPU-parallelizable Sinkhorn for eASOT problem

We have already introduced the definition of the ASOT problem on a single pair
of distributions in the previous subsection. In the case of GPU
parallelization on a multiple OT problem set, we consider computing their
$1$-Wasserstein distances by solving the translated eASOT problems over a
shared anchor space of $\mathcal{W}$ and metric function $d_{\rm AS}$. In this
case, all of the distribution pairs in the batch share a set of parameters of
$\mathcal{W},P$ and $d_{\rm AS}$, which only need to be learned once.
According to Remark 4.16 in [2], we use the GPU parallelization version of the
Sinkhorn’s algorithm to solve the eASOT problem. More concretely, given $N$
translated probability distribution pairs of ASOT
$\\{(\mbox{\boldmath$a$}_{i},\mbox{\boldmath$b$}_{i})\\}_{i=1}^{N}$, where
$\mbox{\boldmath$a$}_{i},\mbox{\boldmath$b$}_{i}\in\mathbb{R}^{k}$, and
$\mbox{\bf
A}:=[\mbox{\boldmath$a$}_{1},\cdots,\mbox{\boldmath$a$}_{N}]\in\mathbb{R}^{k\times
N}$, $\mbox{\bf
B}:=[\mbox{\boldmath$b$}_{1},\cdots,\mbox{\boldmath$b$}_{N}]\in\mathbb{R}^{k\times
N}$, the following equation updates the dual variables
$\mbox{\boldmath$u$}_{i}\in\mathbb{R}^{k}$ and
$\mbox{\boldmath$v$}_{i}\in\mathbb{R}^{k}$ at the $j$-th iteration of
Sinkhorn’s algorithm with GPU parallelization.

$\displaystyle\mbox{\bf U}^{(j+1)}=\frac{\mbox{\bf A}}{\mbox{\bf K}\mbox{\bf
V}^{(j)}},\quad\mbox{\bf V}^{(j+1)}=\frac{\mbox{\bf B}}{\mbox{\bf
K}^{T}\mbox{\bf U}^{(j+1)}},$

where $\mbox{\bf
U}:=[\mbox{\boldmath$u$}_{1},.\cdots,\mbox{\boldmath$u$}_{N}],\mbox{\bf
V}:=[\mbox{\boldmath$v$}_{1},.\cdots,\mbox{\boldmath$v$}_{N}]\in\mathbb{R}^{k\times
N}$, and $\mbox{\bf K}:=\exp\left(-\frac{{\bf
C}_{s}}{\varepsilon}\right)\in\mathbb{R}^{k\times k}$ is the Gibbs kernel.

###  3.4 Complexity analysis

We present a comparison of the time and space complexity of the one-by-one
computation and the batch processing of our ASOT problem with the Sinkhorn’s
solver. Let $N$ denote the number of distributions in the mini-batch, and let
$s$ denote the average number of samples in a distribution.

Time complexity. For the preparation procedure before Sinkhorn’s algorithm,
the one-by-one computation requires an $O(Ns^{2})$ computation to instantiate
kernels, while ASOT only requires $O(Ns)$ to map distributions with mapping
function $P$. For the Sinkhorn’s procedure, the one-by-one computation
requires an $O(s^{3})$ to compute K$v$ and $\mbox{\bf
K}^{T}\mbox{\boldmath$u$}$ for each single Sinkhorn loop. ASOT, on the other
hand, requires $O(k^{3})$, where $k$ denotes the number of anchor points. It
should be noted that the Sinkhorn’s algorithm with block diagonal stacking has
an extra operation to extract the results from the block diagonal matrix with
batch indices in the case of different distribution sizes. On the other hand,
our proposed ASOT with the stacking methods described in Section 3.3 does not
require such operation.

Space complexity. The one-by-one computation requires an $O(Ns^{2})$ space for
restoring every cost matrix, while ASOT only requires $O(k^{2})$. Therefore,
with setting $k$ and $s$ at the same level of magnitude, the ASOT will save
much space and make it able to create a mini-batch of a larger size, which can
also accelerate computation under a limited memory size.

##  4 Anchor Space Learning

As mentioned above, to translate OT problems into ASOT problems, we need to
learn an anchor space to minimize the upper bound in Proposition 1 w.r.t.
$\mathcal{W},P$ and $d_{\rm AS}$. This section presents a deep metric learning
framework, a $k$-means-based framework, and a deep dictionary learning
framework to solve this optimization problem. Each has different advantages
and disadvantages, which will be discussed in this section.

###  4.1 Deep metric learning ASOT (ASOT-ML)

For the framework of our metric learning ASOT, to minimize
$\|\widehat{\mbox{\bf C}}-\mbox{\bf C}\|_{1}$, we first assume that
$\widehat{\mbox{\bf C}}$ is computed by a parameterized metric function
$f(\mathcal{W},P,d_{\rm AS})$, defined as follows.

$\displaystyle\widehat{\mbox{\bf C}}(i,j)$ $\displaystyle=$ $\displaystyle
f(\mathcal{W},P,d_{\rm AS})$ $\displaystyle=$
$\displaystyle\sum_{u,v\in[k]}\mbox{\bf Z}_{x}(i,u)\cdot\mbox{\bf
Z}_{y}(j,u)\cdot d_{\rm AS}(\mbox{\boldmath$w$}_{u},\mbox{\boldmath$w$}_{v})$
$\displaystyle=$
$\displaystyle\sum_{u,v\in[k]}P(\mbox{\boldmath$x$}_{i})(u)\cdot
P(\mbox{\boldmath$y$}_{j})(v)\cdot d_{\rm
AS}(\mbox{\boldmath$w$}_{u},\mbox{\boldmath$w$}_{v}).$

Then the problem becomes a metric learning problem from $f$ to $d_{\rm S}$,
which is the ground metric function of C. Because we have two variables $P$
and $d_{\rm AS}$ that are functions, we parameterize $P$ with a multilayer
perception network (MLP); we use the Mahalanobis distance for $d_{\rm AS}$,
which is computed as $d_{\rm
AS}(\mbox{\boldmath$w$}_{u},\mbox{\boldmath$w$}_{v})=\|\mbox{\bf
M}\mbox{\boldmath$w$}_{u}-\mbox{\bf M}\mbox{\boldmath$w$}_{v}\|_{2}$, where
$\mbox{\bf M}\in\mathbb{R}^{h\times k}$ is a learned parameter matrix, $h$
denotes the number of the hidden dimensions. It should be noted that, in order
to ensure that the outputs of $P$ are within the probability simplex, the
outputs of the MLP for $P$ are normalized by dividing their $l_{1}$ norms.
Thereafter, we adopt an automatic differentiation neural network system to
build our framework, as shown at (a) in Figure 2. The loss function for the
framework is defined as $l(i,j)=|f(\mathcal{W},P,d_{\rm AS})-d_{\rm
S}(\mbox{\boldmath$x$}_{i},\mbox{\boldmath$y$}_{j})|$ for each pair of
sampling points. To avoid the subgradient, we adopt an $l_{2}$ loss as $l_{\rm
ML}=\frac{C}{nm}\sum_{i\in[n],j\in[m]}l(i,j)^{2}$, where $C$ is a scalar.

However, there exist several limitations for ASOT-ML: (i) Quadratic
complexity: Because ASOT-ML needs pairwise distance errors, it has both time
complexity and space complexity of $O(n^{2})$. Therefore, it is impossible to
conduct full pairwise distance computation on a large-scale dataset but have
to use mini-batch pairwise distances, which might affect its performance. (ii)
Unsafe solution: ASOT-ML can not control the scale of the learned cost matrix,
which might cause a risk of overflow while computing Sinkhorn’s algorithm. It
is because when $\mbox{\bf C}_{s}$ has a too large or too small value, the
computation of $\mbox{\bf K}=\exp(-\frac{\mbox{\bf C}_{s}}{\varepsilon})$ will
produce outliers and thus affects the updating process.

###  4.2 $k$-means-based (ASOT-$k$) and Deep dictionary learning ASOT (ASOT-
DL)

Upper bound under the one-hot condition. In order to avoid the limitations of
ASOT-ML, we propose a relative lightweight deep dictionary learning framework
based on several special conditions, which only has a complexity of $O(n)$
level. We present the following proposition of upper bound and provide its
proof sketch, while the complete proof is in the supplementary file.

######  Proposition 2 (Upper bound of absolute 1-Wasserstein error under the
one-hot condition).

Let $\mbox{\boldmath$q$}_{a}^{b}\in\mathbb{R}^{a}$ denotes an $a$-dimensional
one-hot vector of which the $b$-th element is equal to $1$, i.e.
$\mbox{\boldmath$q$}_{a}^{b}(b)=1$. Let
$P_{o}:\mathbb{R}^{d}\to\\{\mbox{\boldmath$q$}_{k}^{i}\\}_{i=1}^{k}$ be an
one-hot mapping function, let
$\mathcal{W}=\\{\mbox{\boldmath$w$}_{i}\\}_{i=1}^{k}$ and $\mbox{\bf
W}=[\mbox{\boldmath$w$}_{1},\cdots,\mbox{\boldmath$w$}_{k}]\in\mathbb{R}^{d\times
k}$. Given two discrete distributions
$\mathcal{D}_{x}=(\mathcal{X},\mbox{\boldmath$a$}),\mathcal{D}_{y}=(\mathcal{Y},\mbox{\boldmath$b$})$,
where
$\mathcal{X}=\\{\mbox{\boldmath$x$}_{i}\\}_{i=1}^{n},\mathcal{Y}=\\{\mbox{\boldmath$y$}_{j}\\}_{j=1}^{m}$,
$\mbox{\boldmath$x$}_{i},\mbox{\boldmath$y$}_{j}\in\mathbb{R}^{d}$, and
$\mbox{\boldmath$a$}\in\Delta^{p}_{n},\mbox{\boldmath$b$}\in\Delta^{p}_{m}$.
Assume an original OT problem $W_{1}(\mathcal{D}_{x},\mathcal{D}_{y},d_{\rm
E})$, which is translated into an ASOT problem as $W_{\rm
AS}(\mathcal{D}_{x},\mathcal{D}_{y},\mathcal{W},P_{o},d_{\rm E})$, where
$d_{\rm E}$ is the Euclidean distance function. Then the following inequality
holds.

$\displaystyle|W_{\rm
AS}(\mathcal{D}_{x},\mathcal{D}_{y},\mathcal{W},P_{o},d_{\rm
E})-W_{1}(\mathcal{D}_{x},\mathcal{D}_{y},d_{\rm
E})|<m\sum_{i\in[n]}\|\mbox{\boldmath$\epsilon$}_{i}^{x}\|_{2}+n\sum_{j\in[m]}\|\mbox{\boldmath$\epsilon$}_{j}^{y}\|_{2},$

where $\mbox{\boldmath$\epsilon$}_{i}^{x}=\mbox{\bf
W}P_{o}(\mbox{\boldmath$x$}_{i})-\mbox{\boldmath$x$}_{i},\mbox{\boldmath$\epsilon$}_{j}^{y}=\mbox{\bf
W}P_{o}(\mbox{\boldmath$y$}_{j})-\mbox{\boldmath$y$}_{j}$ denote the offsets
of reconstructed features to the original ones.

###### Proof.

Let $\eta:=|W_{\rm AS}-W_{1}|$, where $W_{\rm AS}:=W_{\rm
AS}(\mathcal{D}_{x},\mathcal{D}_{y},\mathcal{W},P_{o},d_{\rm E})$, and
$W_{1}:=W_{1}(\mathcal{D}_{x},\mathcal{D}_{y},d_{\rm E})$ . From Proposition
1, we have

$\displaystyle\eta$ $\displaystyle<$ $\displaystyle\|\widehat{\mbox{\bf
C}}-\mbox{\bf C}\|_{1}$ $\displaystyle<$
$\displaystyle\sum_{i\in[n],j\in[m]}\left|\widehat{\mbox{\bf
C}}(i,j)-\mbox{\bf C}(i,j)\right|$ $\displaystyle<$
$\displaystyle\sum_{i\in[n],j\in[m]}\left|\sum_{u,v\in[k]}\mbox{\boldmath$z$}^{x}_{i}(u)\cdot\mbox{\boldmath$z$}^{y}_{j}(v)\cdot\mbox{\bf
C}_{s}(u,v)-\mbox{\bf C}(i,j)\right|,$

where
$\mbox{\boldmath$z$}^{x}_{i}:=P_{o}(\mbox{\boldmath$x$}_{i}),\mbox{\boldmath$z$}^{y}_{j}:=P_{o}(\mbox{\boldmath$y$}_{j})$.
Because $P_{o}$ is an one-hot mapping function, we can let the index of non-
zero element of $P_{o}(\mbox{\boldmath$x$}_{i})$ be $u_{i}$, and let the one
of $P_{o}(\mbox{\boldmath$x$}_{i})$ be $v_{i}$. Substituting this into the
inequalities above obtains

$\eta<\sum_{i\in[n],j\in[m]}\left|\mbox{\bf C}_{s}(u_{i},v_{j})-\mbox{\bf
C}(i,j)\right|.$

From Definition 3, we have $\mbox{\bf C}_{s}(u_{i},v_{j})=d_{\rm
E}(\mbox{\boldmath$w$}_{u_{i}},\mbox{\boldmath$w$}_{v_{j}})=\|\mbox{\boldmath$w$}_{u_{i}}-\mbox{\boldmath$w$}_{v_{j}}\|_{2}$
and $\mbox{\bf C}(i,j)=d_{\rm
E}(\mbox{\boldmath$x$}_{i},\mbox{\boldmath$y$}_{j})=\|\mbox{\boldmath$x$}_{i}-\mbox{\boldmath$y$}_{j}\|_{2}$.
Let $\epsilon_{i}^{x}:=\mbox{\bf
W}P_{o}(\mbox{\boldmath$x$}_{i})-\mbox{\boldmath$x$}_{i},\epsilon_{j}^{y}:=\mbox{\bf
W}P_{o}(\mbox{\boldmath$y$}_{j})-\mbox{\boldmath$y$}_{j}$. Then,
$\epsilon_{i}^{x}=\mbox{\boldmath$w$}_{u_{i}}-\mbox{\boldmath$x$}_{i},\epsilon_{j}^{y}=\mbox{\boldmath$w$}_{v_{j}}-\mbox{\boldmath$y$}_{j}$.
Substituting these into the inequality above obtains

$\eta<\sum_{i\in[n],j\in[m]}\left|\|\mbox{\boldmath$w$}_{u_{i}}-\mbox{\boldmath$w$}_{v_{j}}\|_{2}-\|(\mbox{\boldmath$w$}_{u_{i}}-\mbox{\boldmath$w$}_{v_{j}})+(\mbox{\boldmath$\epsilon$}_{j}^{y}-\mbox{\boldmath$\epsilon$}_{i}^{x})\|_{2}\right|.$

According to the reverse triangle inequality [22], we have

$\displaystyle\eta$ $\displaystyle<$
$\displaystyle\sum_{i\in[n],j\in[m]}\|\mbox{\boldmath$\epsilon$}_{j}^{y}-\mbox{\boldmath$\epsilon$}_{i}^{x}\|_{2}$
$\displaystyle<$
$\displaystyle\sum_{i\in[n],j\in[m]}\|\mbox{\boldmath$\epsilon$}_{j}^{y}\|_{2}+\|\mbox{\boldmath$\epsilon$}_{i}^{x}\|_{2}$
$\displaystyle<$ $\displaystyle
m\sum_{i\in[n]}\|\mbox{\boldmath$\epsilon$}_{i}^{x}\|_{2}+n\sum_{j\in[m]}\|\mbox{\boldmath$\epsilon$}_{j}^{y}\|_{2},$

∎

From Proposition 2, with a special one-hot encoding $P_{o}$ and a shared
$l_{p}$ norm distance function $d_{\rm S}$ for both anchor space and original
space, the distance error could be upper bounded within a summation of
reconstruction errors. In this case, there are only two parameters to be
learned, which are $\mathcal{W}$ and $P_{o}$. Both our proposed ASOT-$k$ and
ASOT-DL methods are based on the same optimization problem derived from
Proposition 2, which is $\mathop{\rm
minimize}_{\mathcal{W},P_{o}}\sum_{i\in[n]}\|\mbox{\boldmath$\epsilon$}_{i}^{x}\|_{2}+\sum_{j\in[m]}\|\mbox{\boldmath$\epsilon$}_{j}^{y}\|_{2}$.
This optimization problem only needs computation of $O(n+m)$ complexity.
Furthermore, for the limitation of the unsafe solution, because of the one-hot
constraint and the share of the metric function $d_{\rm S}$ over both the
anchor space and the original space, the learned cost matrix of anchor points
$\mbox{\bf C}_{s}$ usually has the same scale as the one of original OT
problem C. Therefore, it is less likely to provide unsafe solutions with
Proposition 2. The last, because Proposition 2 is built upon extra conditions
from Proposition 1, its solution set should be included in the later one.
Therefore, theoretically, ASOT-ML can have better solutions than ASOT-$k$ and
ASOT-DL.

By minimizing the reconstruction errors, the problem becomes a classical
clustering problem or a dictionary learning problem. For a clustering problem,
it could be easily solved with the $k$-means algorithm and its variants, where
$\mathcal{W}$ is the set of cluster centers and $P_{o}$ is the clustering
operator. Therefore, we utilize the $k$-means algorithm to learn the anchor
space in our ASOT-$k$ method. However, the one-hot condition might be too
strict, which makes it less possible to achieve low reconstruction errors.
Therefore, We also propose a deep dictionary learning framework ASOT-DL to
optimize a near-one-hot mapping function. In more detail, we adopt a framework
from [23], named deep K-SVD, with adding some modifications to fit our
conditions.

Total loss for deep dictionary learning. We first elaborate on the definition
of the total loss for our deep dictionary learning framework before the
details of the network structure. Let
$\mbox{\boldmath$x$}_{i}\in\mathbb{R}^{d}$ be a $d$-dimensional sampling
point, $\mbox{\bf
X}:=[\mbox{\boldmath$x$}_{1},\mbox{\boldmath$x$}_{2},\cdots,\mbox{\boldmath$x$}_{n}]^{T}\in\mathbb{R}^{n\times
d}$. Let
$\mbox{\boldmath$z$}_{i}=P_{o}(\mbox{\boldmath$x$}_{i})\in\mathbb{R}^{k}$ be a
near-one-hot encoding for $\mbox{\boldmath$x$}_{i}$, $\mbox{\bf
Z}:=[\mbox{\boldmath$z$}_{1},\cdots,\mbox{\boldmath$z$}_{n}]^{T}\in\mathbb{R}^{n\times
k}$. Given a anchor point matrix $\mbox{\bf
W}=[\mbox{\boldmath$w$}_{1},\cdots,\mbox{\boldmath$w$}_{k}]\in\mathbb{R}^{d\times
k}$, we first define a total loss $l$ for our deep dictionary learning
framework as $l(\mbox{\bf W},\mbox{\bf Z}):=\alpha\cdot l_{\rm rc}+\beta\cdot
l_{\rm lp}+\gamma\cdot l_{\rm sp}$, where $l_{\rm rc}$, $l_{\rm lp}$, and
$l_{\rm sp}$ are a reconstruction loss, an $l_{p}$-ball loss, and a simplex
constraint violation loss, respectively.
$\alpha,\beta,\gamma\in\mathbb{R}_{+}$ denote the weights of their
corresponding losses. These losses are computed as $l_{\rm
rc}:=\frac{1}{2n}\sum_{i=1}^{n}\|\mbox{\bf W}(\mbox{\bf
Z}(i,:))^{T}-(\mbox{\bf X}(i,:))^{T}\|_{2}^{2},l_{\rm
lp}:=\frac{1}{n}\sum_{i=1}^{n}\left(\|\mbox{\bf
Z}(i,:)\|_{p}-1\right)^{2},l_{\rm
sp}:=\frac{1}{n}\sum_{i=1}^{n}\left(\mbox{\bf
Z}(i,:)\mbox{\boldmath$1$}_{k}-1\right)^{2},$ where $p>1$ is a scalar. For the
$l_{p}$-ball loss and the simplex constraint violation loss, they penalize the
learned near-one-hot codings to constrain them within an intersection of the
surface of a unit $l_{p}$-ball ($p>1$) and a unit simplex $\Delta^{u}_{k}$,
which is the one-hot vector set $\\{\mbox{\boldmath$q$}_{k}^{i}\\}_{i=1}^{k}$.

Details of the framework structure. In the original paper of deep K-SVD, they
proposed a deep dictionary learning model based on the ISTA algorithm [24]
with an adapted $\lambda$, which is actually a sparse encoder. The ISTA
algorithm solves the lasso problem $\mathop{\rm
minimize}_{\bm{z}}\frac{1}{2}\|\mbox{\bf
W}\mbox{\boldmath$z$}-\mbox{\boldmath$x$}\|_{2}^{2}+\lambda\|\mbox{\boldmath$z$}\|_{1}$
with the proximal gradient descent method that projects the solution of each
iteration into an $l_{1}$-ball by adopting the soft threshold function.
However, because our Definition 3 requires a simplex mapping function $P$, we
shrink the solution domain from an $l_{1}$-ball to a unit simplex
$\Delta^{u}_{k}$ to remove unexpected solutions for faster computation. This
derives the following problem

$\mathop{\rm
minimize}_{{\bm{z}\in\mathbb{R}^{k}_{+}}}~{}\frac{1}{2}\|\mbox{\bf
W}\mbox{\boldmath$z$}-\mbox{\boldmath$x$}\|_{2}^{2}+\lambda(\mbox{\boldmath$z$}^{T}\mbox{\boldmath$1$}_{k}-1),$
(9)

where $\lambda\in\mathbb{R}_{+}$. This problem can be solved by the projected
gradient descent with the orthogonal projection operator as

$\mbox{\boldmath$z$}\leftarrow S_{\lambda}\left(\mbox{\boldmath$z$}-\mbox{\bf
W}^{T}(\mbox{\bf W}\mbox{\boldmath$z$}-\mbox{\boldmath$x$})\right).$ (10)

where $S_{\theta}(\mbox{\boldmath$a$}):\mathbb{R}^{k}\to\mathbb{R}_{+}^{k}$
for a vector $\mbox{\boldmath$a$}\in\mathbb{R}^{k}$ is the orthogonal
projection operator onto the unit simplex, which is defined as
$S_{\theta}(\mbox{\boldmath$a$}):={\rm
max}(\mbox{\boldmath$0$},\mbox{\boldmath$a$}-\theta)$. We present the
derivation of this operator in the supplementary file. Figure 2 (b) shows the
overview of our framework. In our framework, Eq. (10) is represented as a
near-one-hot encoding layer. Our near-one-hot encoder for $P_{o}$ consists of
$L$ sequentially connected near-one-hot encoding layers to generate $z$, where
$L$ denotes the number of layers. We also include a multilayer perception
network (MLP) to evaluate $\lambda$ for each input as deep K-SVD does.

Figure 2: The overall framework of our proposed deep anchor space learning
methods, where $\mbox{\bf X}\in\mathbb{R}^{n\times d}$ denotes the input
feature matrix of sampling points, $\mbox{\bf Z}=P(\mbox{\bf X})$ denotes the
coding matrix mapped by $P$, $\mbox{\bf W}\in\mathbb{R}^{k\times d}$ denotes
the anchor point matrix and $\mbox{\bf M}\in\mathbb{R}^{h\times k}$ denotes
the parameter matrix for metric learning. The mapping function $P$ of ASOT-ML
from (a) is parameterized with an MLP with output normalization, while $P$ of
ASOT-DL from (b) consists of several near-one-hot encoding layers, each of
them computes one iteration of Eq. (10) with an adapted $\lambda$ evaluated by
an MLP.

##  5 Numerical Evaluation

We conduct numerical experiments through three aspects: (i) distance
approximation evaluation, (ii) computational time evaluation, and (iii)
ablation studies of main parameter $k$. This section presents the our
experimental reports.

Figure 3: Average approximation error curves of graph datasets. Table 1:
Average approximation errors of graph datasets (RMSE).

Methods | MUTAG | PTC-MR | ENZYMES | BZR | COX2 | NCI1 | ZINC | MNIST
---|---|---|---|---|---|---|---|---
SW (CPU) | 0.4826 | 1.0729 | 0.3253 | 0.3305 | 0.3253 | 0.6061 | 0.5761 | 0.1084
eOT (CPU) | 0.1155 | 0.2489 | 0.1556 | 0.1366 | 0.0827 | 0.0791 | 0.0754 | 0.0981
BDS-eOT (GPU) | 0.1035 | 0.0675 | 0.1487 | 0.1359 | 0.0807 | 0.0771 | 0.0744 | -
ASOT-ML (CPU) | 0.0264 | 0.0736 | 0.0554 | 0.0144 | 0.0058 | 0.0168 | 0.0306 | 0.0860
ASOT-DL (CPU) | 0.0292 | 0.1037 | 0.0724 | 0.0296 | 0.0293 | 0.0293 | 0.0747 | 0.0972
ASOT-$k$ (CPU) | 0.0862 | 0.1497 | 0.0287 | 0.0285 | 0.0256 | 0.0808 | 0.0410 | 0.0112
eASOT-ML (GPU) | 0.1213 | 0.2211 | 0.1793 | 0.1367 | 0.0803 | 0.0810 | 0.0871 | 0.0766
eASOT-DL (GPU) | 0.1160 | 0.2445 | 0.1685 | 0.1333 | 0.0905 | 0.0965 | 0.0773 | 0.0990
eASOT-$k$ (GPU) | 0.1591 | 0.2318 | 0.1752 | 0.1476 | 0.0940 | 0.1229 | 0.0952 | 0.0766

Table 2: Average computational time cost of graph datasets (sec). For ASOT
methods, the value at the left of “+” denotes the average training time, the
one on the right denotes the distance computation time.

Methods | MUTAG | PTC-MR | ENZYMES | BZR | COX2 | NCI1 | ZINC | MNIST
---|---|---|---|---|---|---|---|---
OT-EMD (CPU) | 2.00 | 6.72 | 27.54 | 16.36 | 28.12 | 1717 | 2566 | 190.4
SW (CPU) | 2.59 | 13.45 | 23.53 | 91.45 | 62.05 | 4512 | 5190 | 25.92
eOT (CPU) | 40.73 | 123.1 | 392.2 | 173.1 | 242.1 | 18574 | 29090 | 572.5
BDS-eOT (GPU) | 6.78 | 26.34 | 85.36 | 42.48 | 56.42 | 3970 | 5703 | -
ASOT-ML (CPU) | 1.98+2.28 | 2.59+18.67 | 1.52+191.4 | 8.01+21.87 | 12.83+27.86 | 3.54+6871 | 1.95+2079 | 0.35+456.8
ASOT-DL (CPU) | 5.97+1.84 | 6.07+7.20 | 5.86+24.46 | 5.12+10.77 | 4.78+21.34 | 7.09+1259 | 0.41+1458 | 0.96+366.1
ASOT-$k$ (CPU) | 0.08+2.66 | 0.10+8.28 | 0.18+32.28 | 0.20+12.98 | 0.21+12.65 | 0.87+1103 | 33.22+1435 | 0.008+88.65
eASOT-ML (GPU) | 1.98+0.03 | 2.59+0.31 | 1.52+2.46 | 8.01+0.20 | 12.83+0.28 | 3.54+93.63 | 1.95+16.43 | 0.35+9.36
eASOT-DL (GPU) | 5.97+0.08 | 6.07+0.38 | 5.86+2.23 | 5.12+0.50 | 4.78+0.32 | 7.09+18.30 | 0.41+16.60 | 0.96+9.31
eASOT-$k$ (GPU) | 0.08+0.64 | 0.10+0.88 | 0.18+3.15 | 0.20+0.91 | 0.21+0.89 | 0.87+86.44 | 33.22+16.96 | 0.008+9.91

###  5.1 Experimental setups

Datasets. We report our experiments on seven real-world graph datasets: MUTAG
[25], PTC-MR [26], ENZYMES [27], BZR, COX2 [28], NCI1[29] and ZINC [30] from
the TUdataset [31]. ZINC is a large-scale dataset containing $249456$ graphs.
For detailed statistical information about each graph dataset, please refer to
the supplementary file. Futhermore, we also report experiments on a computer
vision dataset: MNIST [32]. The MNIST database contains two splits: train and
test subsets, where the train subset has $60000$ images of handwritten digits
from $0$ to $9$, while the test subset has $10000$ images. Each of these
images is resized into a $28\times 28$ matrix.

Preprocessing on graph datasets. Because the node features of the datasets are
one-hot vectors of discrete labels, which are too simple for evaluating
approximation performance. Therefore, to generate more complicated features,
we conduct a feature updating preprocessing before training. More concretely,
We conduct a classic message–passing–based updating equation as follows to
update node features:
$\mbox{\boldmath$x$}_{i}^{\prime}=(1+\epsilon)\cdot\mbox{\boldmath$x$}_{i}+\sum_{j\in\mathcal{N}(i)}\mbox{\boldmath$x$}_{j}$,
where $\mathcal{N}(i)$ denotes the set of indices of adjacent node of $i$-th
node, $\mbox{\boldmath$x$}_{i}\in\mathbb{R}^{d}$ denotes the feature of $i$-th
node and $\mbox{\boldmath$x$}_{i}^{\prime}\in\mathbb{R}^{d}$ denotes its
updated feature. This is exactly the updating equation of the famous graph
isomorphism network [33] without the MLP layer. Specifically, we set the
$\epsilon=0$. For every data, we compute four times of updating equations
iteratively. The final feature for each node will be the concatenation of its
original feature and outputs of each iteration as
$[\mbox{\boldmath$x$}_{i},\mbox{\boldmath$x$}_{i}^{\prime},\mbox{\boldmath$x$}_{i}^{\prime\prime},\mbox{\boldmath$x$}_{i}^{\prime\prime\prime},\mbox{\boldmath$x$}_{i}^{\prime\prime\prime\prime}]\in\mathbb{R}^{5d}$.

Preprocessing on MNIST dataset. We conduct a preprocessing operation on MNIST
dataset to convert the image data into size-variable distributions. In the
preprocessing, we only sample the pixels with non-zero values from the image.
More concretely, assume $i,j\in[28]$, let $p(i,j)\in[256]$ be the pixel value
of $i$-th row and $j$-th column. Each image data is convert to a discrete
distribution with 3-dimensional sampling points as
$\\{[\frac{i}{28},\frac{j}{28},\frac{p(i,j)}{256}]\in\mathbb{R}^{3}:i,j\in[28],p(i,j)>0\\}$.
Because we only sample pixels with values larger than $0$, sizes of these
distributions are variable.

Baselines and other setups. We set the Earth Movers Distance (OT-EMD) solved
by [34] as the baseline distance, which computes an exact solution of the OT
problem in Eq. (1) in the form of linear programming problem. For our proposed
methods, except for ASOT-ML and ASOT-DL, we also implement a variant whose
anchor space is learned by $k$-means according to Proposition 2. eASOT-ML,
eASOT-DL, and eASOT-k denote the entropic ASOT solved by Sinkhorn’s solver
with GPU-parallelization for ASOT-ML, ASOT-DL and the $k$-means variant,
respectively. ASOT-ML, ASOT-DL, and ASOT-k denote the ones with EMD solver.
For other compared approximate methods, we choose the followings: the entropic
OT with the Sinkhorn’s algorithm [2] on CPU (eOT), the entropic OT with the
Sinkhorn’s algorithm and the block diagonal stacking on GPU (BDS-eOT), and the
sliced Wasserstein distance (SW) [9]. For the parameter settings of other
methods, we use a Sinkhorn’s solver with $\varepsilon=0.1$ and $50$
iterations, a sliced Wasserstein distance solver with $d$ projections. For our
methods, we set the number of anchor points $k$ equal to the maximum number of
nodes that a graph contains of the dataset. For the weights of losses of ASOT-
DL, we set $\alpha=1,\beta=1,\gamma=1$. For ASOT-ML, we set $C=100$. More
detailed parameter settings are provided in the supplementary file. We conduct
$500$ epochs of training with a batch size of $500$ (graphs), except for ZINC,
we only compute 50 epochs. For the subset division of graph datasets, we run a
$10$-fold cross-validation where train vs. test is $9:1$, except for ZINC,
where we use the default splits. We compute the pairwise distance of the whole
dataset (train + test) for evaluation, except for ZINC, where we only compute
the test subset. For the subset division of the MNIST dataset, we randomly
choose $50$ images for each category of digits from the test subset, with
which we construct a subset of $500$ images. We conduct $10$ times of distance
approximation experiments on these subsets, which are randomly generated with
different random seeds for each run. The machine for experiments is equipped
with an Intel CPU i7-12700KF, 32 GB RAM, and a GPU of NVIDIA GeForce RTX
3090Ti.

###  5.2 Wasserstein distance approximation

Distance approximation evaluation. In this experiments, we simulate a general
task of pairwise Wasserstein distance matrix computation, where the distance
between two graphs (images) is defined as the $1$-Wasserstein distance with a
Euclidean ground metric between their node feature sets (pixel feature sets).
For the evaluation criterion, we adopt the rooted mean square error (RMSE)
between approximated distance matrices and the ground truth computed by EMD.
Figures 3 shows the curves of approximation errors of three graph datasets,
where we can see both the eASOT-DL and eASOT-ML show great performance, which
is very close to eOT and BDS-eOT. Table 2 shows the average approximation
errors. We omit the standard deviations here but provide them in the
supplementary file because all of them are less than $0.05$.

Computational time cost evaluation. We also record and evaluate the distance
computation time cost in the distance approximation experiments. Table 2 shows
the time records of the computation of distance matrices in seconds, where the
fastest one is with bold typeface. Our proposed eASOT-ML, eASOT-DL, and
eASOT-k with GPU parallelization show amazing advantages on time cost compared
to other conventional methods. We can observe that, on some datasets, the time
cost of eASOT-ML or ASOT-ML is larger than the one of eASOT-DL or ASOT-DL,
which fits our prediction that ASOT-ML may compute unsafe solutions that
increase the number of required solving iterations. We provide the model
training time costs in the supplementary file, which are also low costs with
GPU parallelization.

Figure 4: Average approximation error and computational time cost curves of an
ablation study of the parameter $k$, where we report the results of a series
of $k$ of $[5,10,20,30,40,50,60,70,80,90,100,200]$.

###  5.3 Ablation studies

As described in previous sections, the parameter $k$ is considered the most
important parameter of our proposed ASOT, which denotes the number of
learnable anchor points. In order to study how the approximation error
deteriorates and computational time decreases with different $k$, we report
two ablation study experiments on two different datasets: the graph dataset
MUTAG and the computer vision dataset MNIST. We prepare a candidate set of $k$
as $[5,10,20,30,40,50,60,70,80,90,100,200]$. For each dataset, we conduct 5
runs of a single experiment and then report the means and standard deviations
of the approxiamtion error. We keep the experimental setups as the same as
described in Section 5.1, except for the number of the Sinkhorn iterations,
which we set as $200$ for MUTAG and $500$ for MNIST. Figure 4 shows the
experimental results of approximation error and computational time with
different $k$. From the curves one can observe a relatively stable descend of
the approximation error curves of eASOT-ML on the both datasets, while the one
of eASOT-DL shows signs of rebound on MUTAG dataset when $k$ is larger than
90. We belive that this is because eASOT-DL optimize a loss which is a upper
bound of the one of eASOT-ML, which make it easier to reach the limits and
thus easier to overfit.

##  6 Conclusion

This paper presents a translation of the OT problem into an ASOT problem to
accelerate batch processing of multiple OT problems, with taking advantages of
the potential common characteristics of distributions. In the proposed ASOT,
we only transport masses over a learned and fixed anchor space. We then prove
the upper bound of the Wasserstein distance approximation error of such ASOT
problem. Based on this, we propose a metric learning (ASOT-ML), a $k$-means-
based (ASOT-$k$) and a deep dictionary learning framework (ASOT-DL) for anchor
space learning. Experimental results demonstrate that our methods
significantly decrease computational time while maintaining a reasonable level
of approximation performance. For future work, we aim to design closed-form
expressions for the mapping function and metric function of ASOT, eliminating
the requirement for anchor space learning.

##  Appendix A Proofs and derivation

This section provides the proofs of Lemma 1, Proposition 1, Proposition 2, and
a derivation of Eq. (9) of the main paper. For mathematical notations, we
follow what we have described in Section 2 of the main paper.

###  A.1 Proof of Lemma 1

###### Proof.

Let $\mbox{\bf C}_{s}\in\mathbb{R}^{k\times k}$ denote a $k\times k$ cost
matrix. Let $\mbox{\bf Z}_{x}\in\mathbb{R}^{n\times k},\mbox{\bf
Z}_{y}\in\mathbb{R}^{m\times k}$, where each row of $\mbox{\bf Z}_{x}$ and
$\mbox{\bf Z}_{y}$ is a simplex vector. Assume an OT problem over discrete
distributions with propbability vectors
$\mbox{\boldmath$a$}\in\Delta^{p}_{n},\mbox{\boldmath$b$}\in\Delta^{p}_{m}$ as
follows.

$W_{1}=\mathop{\rm minimize}_{{\bf
T}\in\mathcal{U}(\bm{a},\bm{b})}\langle\mbox{\bf T},\widehat{\mbox{\bf
C}}\rangle,$

where $\widehat{\mbox{\bf C}}\in\mathbb{R}^{n\times m}$ can be decomposed into
$\widehat{\mbox{\bf C}}=\mbox{\bf Z}_{x}\mbox{\bf C}_{s}\mbox{\bf Z}^{T}_{y}$.
Then we have

$\displaystyle W_{1}$ $\displaystyle=$ $\displaystyle\mathop{\rm
minimize}_{{\bf
T}\in\mathcal{U}(\bf{a},\bf{b})}\sum_{i\in[n],j\in[m]}\mbox{\bf
T}(i,j)\cdot\widehat{\mbox{\bf C}}(i,j)$ $\displaystyle=$
$\displaystyle\mathop{\rm minimize}_{{\bf
T}\in\mathcal{U}(\bf{a},\bf{b})}\sum_{i\in[n],j\in[m]}\mbox{\bf
T}(i,j)\cdot\left(\sum_{u,v\in[k]}\mbox{\bf Z}_{x}(i,u)\mbox{\bf
Z}_{y}(j,v)\mbox{\bf C}_{s}(u,v)\right)$ $\displaystyle=$
$\displaystyle\mathop{\rm minimize}_{{\bf
T}\in\mathcal{U}(\bf{a},\bf{b})}\sum_{i\in[n],j\in[m]}\sum_{u,v\in[k]}\mbox{\bf
Z}_{x}(i,u)\cdot\mbox{\bf Z}_{y}(j,v)\cdot\mbox{\bf T}(i,j)\cdot\mbox{\bf
C}_{s}(u,v)$ $\displaystyle=$ $\displaystyle\mathop{\rm minimize}_{{\bf
T}\in\mathcal{U}(\bf{a},\bf{b})}\sum_{u,v\in[k]}\left(\sum_{i\in[n],j\in[m]}\mbox{\bf
Z}_{x}(i,u)\cdot\mbox{\bf Z}_{y}(j,v)\cdot\mbox{\bf
T}(i,j)\right)\cdot\mbox{\bf C}_{s}(u,v).$

Let $\mbox{\bf T}_{s}=\mbox{\bf Z}^{T}_{x}\mbox{\bf T}\mbox{\bf
Z}_{y}\in\mathbb{R}^{k\times k}$, the above equations can be converted as
follows.

$W_{1}^{u}=\mathop{\rm minimize}_{{\bf
T}\in\mathcal{U}(\bf{a},\bf{b})}\langle\mbox{\bf T}_{s},\mbox{\bf
C}_{s}\rangle.$ (11)

To reparameterize this problem with $\mbox{\bf T}_{s}$, we have to find the
domain of $\mbox{\bf T}_{s}$. Right multiplying $\mbox{\bf T}_{s}$ with
$\mbox{\boldmath$1$}_{k}$ obtains

$\mbox{\bf T}_{s}\mbox{\boldmath$1$}_{k}=\mbox{\bf Z}^{T}_{x}\mbox{\bf
T}\mbox{\bf Z}_{y}\mbox{\boldmath$1$}_{k}.$ (12)

Because each row of $\mbox{\bf Z}_{y}$ is a simplex vector, we have $\mbox{\bf
Z}_{y}\mbox{\boldmath$1$}_{k}=\mbox{\boldmath$1$}_{m}$. Subtituting it into
Eq. (12) obtains

$\mbox{\bf T}_{s}\mbox{\boldmath$1$}_{k}\quad=\quad\mbox{\bf
Z}_{x}^{T}\mbox{\bf T}\mbox{\boldmath$1$}_{m}\quad=\quad\frac{1}{n}\mbox{\bf
Z}^{T}_{x}\mbox{\boldmath$a$}.$

In the same way, we can obtain $\mbox{\bf
T}_{s}^{T}\mbox{\boldmath$1$}_{k}=\frac{1}{m}\mbox{\bf
Z}^{T}_{y}\mbox{\boldmath$b$}$. Futhermore, because for all $i\in[n],j\in[m]$,
we have $\mbox{\bf T}(i,j)\geq 0$, then $\forall u,v\in[k],\mbox{\bf
T}_{s}(u,v)\geq 0$ holds. In summary, $\mbox{\bf
T}_{s}\in\mathcal{U}(\mbox{\boldmath$a$}^{\prime},\mbox{\boldmath$b$}^{\prime})$,
where

$\displaystyle\mbox{\boldmath$a$}^{\prime}=\frac{1}{n}\mbox{\bf
Z}^{T}_{x}\mbox{\boldmath$a$}=\frac{1}{n}\sum_{i=1}^{n}\mbox{\bf
Z}_{x}(i,:)\mbox{\boldmath$a$}(i),\quad\mbox{\boldmath$b$}^{\prime}=\frac{1}{m}\mbox{\bf
Z}^{T}_{y}\mbox{\boldmath$b$}=\frac{1}{n}\sum_{i=j}^{m}\mbox{\bf
Z}_{y}(j,:)\mbox{\boldmath$b$}(j).$

Finally, Eq. (11) is reparameterized as

$W_{1}=\mathop{\rm minimize}_{{\bf
T}_{s}\in\mathcal{U}(\bf{a}^{\prime},\bf{b}^{\prime})}\ \langle\mbox{\bf
T}_{s},\mbox{\bf C}_{s}\rangle.$ (13)

Give $\mathcal{X},\mathcal{Y},\mathcal{W}$, mapping function $P$, and metric
function $d_{\rm AS}$ as defined in Definition 3, which satisfy $\forall
u,v\in[k],\mbox{\bf C}_{s}(u,v)=d_{\rm
AS}(\mbox{\boldmath$w$}_{i},\mbox{\boldmath$w$}_{j})$ and $\forall
i\in[n],\mbox{\bf Z}_{x}(i,:)=P(\mbox{\boldmath$x$}_{i})$ and $\forall
j\in[m],\mbox{\bf Z}_{x}(j,:)=P(\mbox{\boldmath$x$}_{j})$. Then Eq. (13) is an
ASOT problem of $W_{\rm
AS}(\mathcal{D}_{x},\mathcal{D}_{y},\mathcal{W},P,d_{\rm AS})$, which
completes the proof. ∎

###  A.2 Proof of Proposition 1

###### Proof.

According to Lemma 1, the ASOT problem of $W_{\rm
AS}(\mathcal{D}_{x},\mathcal{D}_{y},\mathcal{W},P,d_{\rm AS})$ is equivalant
to the OT problem defined in Eq. (7). Let $\widehat{\mbox{\bf T}}^{\star}$ be
the optimal parameter of the equivalant OT problem in Eq. (7), let $\mbox{\bf
T}^{\star}$ be the optimal parameter of the original OT problem of
$W_{1}(\mathcal{D}_{x},\mathcal{D}_{y},d_{\rm S})$. we have

$\displaystyle|W_{\rm AS}-W_{1}|$ $\displaystyle:=$ $\displaystyle|W_{\rm
AS}(\mathcal{D}_{x},\mathcal{D}_{y},\mathcal{W},P,d_{\rm
AS})-W_{1}(\mathcal{D}_{x},\mathcal{D}_{y},d_{\rm S})|$ $\displaystyle=$
$\displaystyle\left|\langle\widehat{\mbox{\bf T}}^{\star},\widehat{\mbox{\bf
C}}\rangle-\langle\mbox{\bf T}^{\star},\mbox{\bf C}\rangle\right|$
$\displaystyle=$ $\displaystyle\left|\sum_{i\in[n],j\in[m]}\widehat{\mbox{\bf
T}}^{\star}(i,j)\widehat{\mbox{\bf C}}(i,j)-\sum_{i\in[n],j\in[m]}\mbox{\bf
T}^{\star}(i,j)\mbox{\bf C}(i,j)\right|.$

In the case of $W_{\rm AS}\mathcal{D}_{x},\mathcal{D}_{y},\mathcal{W},P,d_{\rm
AS})>W_{1}(\mathcal{D}_{x},\mathcal{D}_{y},d_{\rm S})$, we have

$\displaystyle|W_{\rm AS}-W_{1}|$ $\displaystyle=$
$\displaystyle\sum_{i\in[n],j\in[m]}\widehat{\mbox{\bf
T}}^{\star}(i,j)\widehat{\mbox{\bf C}}(i,j)-\sum_{i\in[n],j\in[m]}\mbox{\bf
T}^{\star}(i,j)\mbox{\bf C}(i,j)$ $\displaystyle\leq$
$\displaystyle\sum_{i\in[n],j\in[m]}\mbox{\bf
T}^{\star}(i,j)\widehat{\mbox{\bf C}}(i,j)-\sum_{i\in[n],j\in[m]}\mbox{\bf
T}^{\star}(i,j)\mbox{\bf C}(i,j)$ $\displaystyle\leq$
$\displaystyle\sum_{i\in[n],j\in[m]}\mbox{\bf
T}^{\star}(i,j)\left(\widehat{\mbox{\bf C}}(i,j)-\mbox{\bf C}(i,j)\right),$

which is because OT problem is convex, its optimal solution is the global
minima. Moreover, because $\mbox{\bf
T}^{\star}\in\mathcal{U}(\mbox{\boldmath$a$},\mbox{\boldmath$b$})$ and
$\mbox{\boldmath$a$},\mbox{\boldmath$b$}$ are both in the probability
simplexes, we have $\sum_{i\in[n],j\in[m]}\mbox{\bf
T}^{\star}(i,j)=\|\mbox{\boldmath$a$}\|_{1}=\|\mbox{\boldmath$b$}\|_{1}=1$.
Then we have

$\displaystyle|W_{\rm AS}-W_{1}|$ $\displaystyle\leq$
$\displaystyle\sum_{i\in[n],j\in[m]}\mbox{\bf
T}^{\star}(i,j)\left(\widehat{\mbox{\bf C}}(i,j)-\mbox{\bf C}(i,j)\right)$
$\displaystyle\leq$ $\displaystyle\sum_{i\in[n],j\in[m]}\mbox{\bf
T}^{\star}(i,j)\left|\widehat{\mbox{\bf C}}(i,j)-\mbox{\bf C}(i,j)\right|$
$\displaystyle<$ $\displaystyle\|\widehat{\mbox{\bf C}}-\mbox{\bf
C}\|_{\infty}$ $\displaystyle<$ $\displaystyle\|\widehat{\mbox{\bf
C}}-\mbox{\bf C}\|_{1}.$

In the case of $W_{\rm
AS}(\mathcal{D}_{x},\mathcal{D}_{y},\mathcal{W},P,d_{\rm AS})\leq
W_{1}(\mathcal{D}_{x},\mathcal{D}_{y},d_{\rm S})$, we have

$\displaystyle|W_{\rm AS}-W_{1}|$ $\displaystyle=$
$\displaystyle\sum_{i\in[n],j\in[m]}\mbox{\bf T}^{\star}(i,j)\mbox{\bf
C}(i,j)-\sum_{i\in[n],j\in[m]}\widehat{\mbox{\bf
T}}^{\star}(i,j)\widehat{\mbox{\bf C}}(i,j)$ $\displaystyle\leq$
$\displaystyle\sum_{i\in[n],j\in[m]}\widehat{\mbox{\bf
T}}^{\star}(i,j)\left(\mbox{\bf C}(i,j)-\widehat{\mbox{\bf C}}(i,j)\right)$
$\displaystyle\leq$ $\displaystyle\sum_{i\in[n],j\in[m]}\widehat{\mbox{\bf
T}}^{\star}(i,j)\left|\mbox{\bf C}(i,j)-\widehat{\mbox{\bf C}}(i,j)\right|$
$\displaystyle<$ $\displaystyle\|\mbox{\bf C}-\widehat{\mbox{\bf
C}}\|_{\infty}$ $\displaystyle<$ $\displaystyle\|\mbox{\bf
C}-\widehat{\mbox{\bf C}}\|_{1}.$

Combining the results of both cases above completes the proof. ∎

###  A.3 Proof of Proposition 2

Before proving Proposition 2, we first prove the following theorem.

######  Theorem 1.

The function of $l_{2}$ norm on an $\mathbb{R}^{d}$ vector space, denoted as
$f(\mbox{\boldmath$x$}):=\|\mbox{\boldmath$x$}\|_{2}$, is $1$-Lipschitz
continuous w.r.t. the Euclidean metric.

###### Proof.

According to the reverse triangle inequality [22] we have

$\displaystyle|f(\mbox{\boldmath$x$})-f(\mbox{\boldmath$y$})|$
$\displaystyle=$
$\displaystyle\left|\|\mbox{\boldmath$x$}\|_{2}-\|\mbox{\boldmath$y$}\|_{2}\right|$
$\displaystyle\leq$
$\displaystyle\|\mbox{\boldmath$x$}-\mbox{\boldmath$y$}\|.$

Therefore, $f(\mbox{\boldmath$x$})$ is $1$-Lipschitz continuous, which
completes the proof. ∎

###### Proof.

According to Proposition 1, we have $|W_{\rm
AS}(\mathcal{D}_{x},\mathcal{D}_{y},\mathcal{W},P_{o},d_{\rm
E})-W_{1}(\mathcal{D}_{x},\mathcal{D}_{y},d_{\rm E})|<\|\widehat{\mbox{\bf
C}}-\mbox{\bf C}\|$. Expanding this inequality obtains

$\displaystyle|W_{\rm AS}-W_{1}|$ $\displaystyle:=$ $\displaystyle|W_{\rm
AS}(\mathcal{D}_{x},\mathcal{D}_{y},\mathcal{W},P_{o},d_{\rm
E})-W_{1}(\mathcal{D}_{x},\mathcal{D}_{y},d_{\rm E})|$ $\displaystyle<$
$\displaystyle\|\widehat{\mbox{\bf C}}-\mbox{\bf C}\|$ $\displaystyle<$
$\displaystyle\sum_{i\in[n],j\in[m]}\left|\widehat{\mbox{\bf
C}}(i,j)-\mbox{\bf C}(i,j)\right|_{1}.$

Because $\widehat{\mbox{\bf C}}$ can be decomposed into $\widehat{\mbox{\bf
C}}=\mbox{\bf Z}_{x}\mbox{\bf C}_{s}\mbox{\bf Z}^{T}_{y}$ according to
Proposition 1, we have

$\displaystyle|W_{\rm AS}-W_{1}|$ $\displaystyle<$
$\displaystyle\sum_{i\in[n],j\in[m]}\left|\widehat{\mbox{\bf
C}}(i,j)-\mbox{\bf C}(i,j)\right|$ $\displaystyle=$
$\displaystyle\sum_{i\in[n],j\in[m]}\left|\sum_{u,v\in[k]}\mbox{\bf
Z}_{x}(i,u)\cdot\mbox{\bf Z}_{y}(j,v)\cdot\mbox{\bf C}_{s}(u,v)-\mbox{\bf
C}(i,j)\right|$ $\displaystyle=$
$\displaystyle\sum_{i\in[n],j\in[m]}\left|\sum_{u,v\in[k]}P_{o}(\mbox{\boldmath$x$}_{i})(u)\cdot
P_{o}(\mbox{\boldmath$y$}_{j})(v)\cdot\mbox{\bf C}_{s}(u,v)-\mbox{\bf
C}(i,j)\right|.$

Because $P_{o}$ is an one-hot mapping function, we can let the index of non-
zero element of $P_{o}(\mbox{\boldmath$x$}_{i})$ be $u_{i}$, and let the one
of $P_{o}(\mbox{\boldmath$x$}_{i})$ be $v_{i}$. Substituting this into the
equation above obtains

$|W_{\rm AS}-W_{1}|<\sum_{i\in[n],j\in[m]}\left|\mbox{\bf
C}_{s}(u_{i},v_{j})-\mbox{\bf C}(i,j)\right|.$

From Definition 3, we have $\mbox{\bf C}_{s}(u_{i},v_{j})=d_{\rm
S}(\mbox{\boldmath$w$}_{u_{i}},\mbox{\boldmath$w$}_{v_{j}})=\|\mbox{\boldmath$w$}_{u_{i}}-\mbox{\boldmath$w$}_{v_{j}}\|_{2}$
and $\mbox{\bf C}(i,j)=d_{\rm
S}(\mbox{\boldmath$x$}_{i},\mbox{\boldmath$y$}_{j})=\|\mbox{\boldmath$x$}_{i}-\mbox{\boldmath$y$}_{j}\|_{2}$.
Let $\epsilon_{i}^{x}:=\mbox{\bf
W}P_{o}(\mbox{\boldmath$x$}_{i})-\mbox{\boldmath$x$}_{i},\epsilon_{j}^{y}:=\mbox{\bf
W}P_{o}(\mbox{\boldmath$y$}_{j})-\mbox{\boldmath$y$}_{j}$. Then,
$\epsilon_{i}^{x}=\mbox{\boldmath$w$}_{u_{i}}-\mbox{\boldmath$x$}_{i},\epsilon_{j}^{y}=\mbox{\boldmath$w$}_{v_{j}}-\mbox{\boldmath$y$}_{j}$.
Substituting these into the inequality above obtains

$\displaystyle|W_{\rm AS}-W_{1}|$ $\displaystyle<$
$\displaystyle\sum_{i\in[n],j\in[m]}\left|\|\mbox{\boldmath$w$}_{u_{i}}-\mbox{\boldmath$w$}_{v_{j}}\|_{2}-\|(\mbox{\boldmath$w$}_{u_{i}}-\mbox{\boldmath$\epsilon$}_{i}^{x})-(\mbox{\boldmath$w$}_{v_{j}}-\mbox{\boldmath$\epsilon$}_{j}^{y})\|_{2}\right|$
$\displaystyle=$
$\displaystyle\sum_{i\in[n],j\in[m]}\left|\|\mbox{\boldmath$w$}_{u_{i}}-\mbox{\boldmath$w$}_{v_{j}}\|_{2}-\|(\mbox{\boldmath$w$}_{u_{i}}-\mbox{\boldmath$w$}_{v_{j}})+(\mbox{\boldmath$\epsilon$}_{j}^{y}-\mbox{\boldmath$\epsilon$}_{i}^{x})\|_{2}\right|.$

According to Theorem 1, we have

$\displaystyle|W_{\rm AS}-W_{1}|$ $\displaystyle<$
$\displaystyle\sum_{i\in[n],j\in[m]}\|\mbox{\boldmath$\epsilon$}_{j}^{y}-\mbox{\boldmath$\epsilon$}_{i}^{x}\|_{2}=\sum_{i\in[n],j\in[m]}\|\mbox{\boldmath$\epsilon$}_{j}^{y}+(-\mbox{\boldmath$\epsilon$}_{i}^{x})\|_{2}$
$\displaystyle\leq$
$\displaystyle\sum_{i\in[n],j\in[m]}\|\mbox{\boldmath$\epsilon$}_{j}^{y}\|_{2}+\|\mbox{\boldmath$\epsilon$}_{i}^{x}\|_{2}=m\sum_{i\in[n]}\|\mbox{\boldmath$\epsilon$}_{i}^{x}\|_{2}+n\sum_{j\in[m]}\|\mbox{\boldmath$\epsilon$}_{j}^{y}\|_{2},$

which completes the proof. ∎

###  A.4 Derivation of Eq. (10)

Recall the primal problem that is expressed as

$\displaystyle\mathop{\rm
minimize}_{\bm{z}\in\mathbb{R}^{k}_{+},\langle\bm{z},\mbox{\boldmath$1$}_{k}\rangle\leq
1}\left\\{f(\mbox{\boldmath$z$}):=\frac{1}{2}\|\mbox{\bf
W}\mbox{\boldmath$z$}-\mbox{\boldmath$x$}\|_{2}^{2}\right\\},$ (14)

where we want to optimize a $z$ on the unit simplex
$\mbox{\boldmath$\Delta$}^{u}_{k}$ that minimizies the squared reconstruction
error. We first ignore the non-negative constraint
$\mbox{\boldmath$z$}\in\mathbb{R}^{k}_{+}$, and we compute the Lagrangian of
$f(\mbox{\boldmath$z$})$ in Eq. (14) as follows.

$\displaystyle L(\mbox{\boldmath$z$},\lambda):=\frac{1}{2}\|\mbox{\bf
W}\mbox{\boldmath$z$}-\mbox{\boldmath$x$}\|_{2}^{2}+\lambda(\langle\mbox{\boldmath$z$},\mbox{\boldmath$1$}_{k}\rangle-1).$

Assuming that the value of $\lambda>0$ is already known, the minimization of
the above Lagrangian with respect to
$\mbox{\boldmath$z$}\in\mathbb{R}^{k}_{+}$ can be expressed as follows.

$\displaystyle\mathop{\rm
minimize}_{\bm{z}\in\mathbb{R}^{k}_{+}}\frac{1}{2}\|\mbox{\bf
W}\mbox{\boldmath$z$}-\mbox{\boldmath$x$}\|_{2}^{2}+\langle\mbox{\boldmath$z$},\lambda\mbox{\boldmath$1$}_{k}\rangle.$
(15)

Let
$g(\mbox{\boldmath$z$})=\langle\mbox{\boldmath$z$},\lambda\mbox{\boldmath$1$}_{k}\rangle$,
according to the proximal gradient descent method [35], the update equation of
$z$ with a fixed step size as $1$ is expressed as follows.

$\displaystyle\mbox{\boldmath$z$}^{(t+1)}$ $\displaystyle=$ $\displaystyle{\rm
prox}_{g}(\mbox{\boldmath$z$}^{(t)}-\nabla f(\mbox{\boldmath$z$}^{(t)}))$ (16)
$\displaystyle=$ $\displaystyle{\rm
prox}_{g}(\mbox{\boldmath$z$}^{(t)}-\mbox{\bf W}^{T}(\mbox{\bf
W}\mbox{\boldmath$z$}^{(t)}-\mbox{\boldmath$x$})),$ (17)

where ${\rm prox}_{g}$ denotes the proximal operator, which is defined as

$\displaystyle{\rm prox}_{g}(\mbox{\boldmath$a$})$
$\displaystyle\\!\\!=\\!\\!$ $\displaystyle\mathop{\rm
arg~{}min}_{\bm{a}}g(\mbox{\boldmath$a$})+\frac{1}{2}\|\mbox{\boldmath$z$}-\mbox{\boldmath$a$}\|_{2}^{2}$
(18) $\displaystyle\\!\\!=\\!\\!$ $\displaystyle\mathop{\rm
arg~{}min}_{\bm{a}}\frac{1}{2}\|\mbox{\boldmath$z$}-\mbox{\boldmath$a$}\|_{2}^{2}+\langle\mbox{\boldmath$a$},\lambda\mbox{\boldmath$1$}_{k}\rangle$
(19) $\displaystyle\\!\\!=\\!\\!$ $\displaystyle\mathop{\rm
arg~{}min}_{\bm{a}}\sum_{i=1}^{k}\left(\frac{1}{2}(\mbox{\boldmath$z$}(i)-\mbox{\boldmath$a$}(i))^{2}+\lambda\mbox{\boldmath$a$}(i)\right).\quad\quad$
(20)

Eq. (20) can be written as $k$ subproblems with respect to each element of
$a$, such as

$\displaystyle{\rm prox}_{g}(\mbox{\boldmath$a$})(i)=\mathop{\rm
arg~{}min}_{\bm{a}(i)}\frac{1}{2}(\mbox{\boldmath$z$}(i)-\mbox{\boldmath$a$}(i))^{2}+\lambda\mbox{\boldmath$a$}(i).$

Because of the non-negative condition of $z$ and Eq. (17), we have $\forall
i\in[k],\mbox{\boldmath$a$}(i)>0$. Therefore, the solutions of the above
subproblems are
$\mbox{\boldmath$a$}(i)=\max(0,\mbox{\boldmath$z$}(i)-\lambda)$. Generalizing
these solutions to the vector $a$ yields

$\displaystyle{\rm
prox}_{g}(\mbox{\boldmath$a$})=\max(\mbox{\boldmath$0$},\mbox{\boldmath$a$}-\lambda\mbox{\boldmath$1$}_{k}),$
(21)

which derives the equation of our orthogonal projection operator
$S_{\theta}(\mbox{\boldmath$a$})$ of the paper. Substituting Eq. (21) into Eq.
(17) yields

$\displaystyle\mbox{\boldmath$z$}^{(t+1)}=S_{\lambda}(\mbox{\boldmath$z$}^{(t)}-\mbox{\bf
W}^{T}(\mbox{\bf W}\mbox{\boldmath$z$}^{(t)}-\mbox{\boldmath$x$})),$ (22)

which completes the derivation.

Table 3: Statistical information of datasets

Dataset | Dim. | Graphs | Avg. $|\mathcal{V}|$ | Avg. $|\mathcal{E}|$ | Max. $|\mathcal{V}|$
---|---|---|---|---|---
MUTAG | 7 | 188 | 17.93 | 19.79 | 28
PTC-MR | 18 | 344 | 14.29 | 14.69 | 64
ENZYMES | 3 | 600 | 32.63 | 62.14 | 126
BZR | 53 | 405 | 35.75 | 38.35 | 57
COX2 | 35 | 467 | 41.22 | 43.45 | 56
NCI1 | 37 | 4110 | 29.86 | 32.30 | 111
ZINC | 28 | 249456 | 23.15 | 24.90 | 38

Table 4: ASOT parameter settings

Parameters | ASOT-ML | ASOT-DL | ASOT-$k$
---|---|---|---
Num. of anchor points $k=|\mathcal{W}|$ | Max. $|\mathcal{V}|$ | Max. $|\mathcal{V}|$ | Max. $|\mathcal{V}|$
Hidden dim. of Mahalanobis distance $h$ | $d$ | - | -
MLP for $P$ of ASOT-ML | $(5d,10d,k)$ | - | -
MLP for $\lambda$ evaluation of ASOT-DL | - | $(5d,2.5d,1)$ | -
Num. of near-one-hot layers | - | 20 | -
$C$ | 100 | - | -
$\alpha$ | - | 1.0 | -
$\beta$ | - | 1.0 | -
$\gamma$ | - | 1.0 | -
Sinkhorn $\varepsilon$ | 0.1 | 0.1 | 0.1
Sinkhorn iterations | 50 | 50 | 50
Optimizer | Adam[36] | Adam[36] | Mini-batch $k$-means [37]
Learning rate | 0.01 | 0.01 | -

##  Appendix B Details of graph experiments

This section describes the detailed settings and results of our numerical
experiments on graph datasets, which are not included in our main paper.

Detailed statistical information of graph datasets. Table 4 records the
statistical information of graph datasets used in our numerical experiments,
where “Dim.” denotes the number of feature dimensions, “Avg. $|\mathcal{V}|$”
and “Avg. $|\mathcal{E}|$” denote the average number of nodes and edges for a
single graph, and “Max. $|\mathcal{V}|$” denotes the maximum number of nodes
for a single graph. Futhermore, ZINC dataset has a default splits of train and
test subsets, whose sizes are 220011 and 5000, respectively.

Detailed ASOT parameter settings. Table 4 records the parameter settings of
ASOT-ML, ASOT-DL, and ASOT-$k$ in our graph experiments, where $d$ denotes the
number of dimensions of node features before preprocessing, Max. $\mathcal{V}$
denotes the same one in Table 4. For the parameters we do not record in the
main paper, we set the hidden dimensions of $\mbox{\bf
M}\in\mathbb{R}^{h\times k}$ equal to $d$ in ASOT-ML. Furthermore, we
parameterize $P$ with a two-layer MLP with entries of $(5d,10d,k)$ in ASOT-ML.
We also use another two-layer MLP with entries of $(5d,2.5d,1)$ to evaluate
$\lambda$ in ASOT-DL.

Standard deviations of approximation errors Table 6 reports the standard
deviations of approximation errors in our graph experiments, w.r.t. the
average errors in Table 1 of the main paper.

Training time costs. Table 6 reports the average training time costs of our
ASOT methods on different graph datasets. Traning time costs in ZINC is less
than other datasets because we only conduct 50 epochs of training on ZINC,
while we train 500 epochs on other datasets. From the table, we can see that
ASOT-ML is faster than ASOT-DL on small datasets such as MUTAG, PTC-MR, and
ENZYMES. This is because ASOT-DL includes $20$ near-one-hot layers for $P$
while ASOT-ML only has a two-layer MLP encoder, which will cost more training
time.

Table 5: Standard deviations of approximation errors of graph datasets (RMSE).

Methods | MUTAG | PTC-MR | ENZYMES | BZR | COX2 | NCI1 | ZINC
---|---|---|---|---|---|---|---
ASOT-ML | 0.0244 | 0.0164 | 0.0555 | 0.0015 | 0.0006 | 0.0022 | 0.0063
ASOT-DL | 0.0104 | 0.0350 | 0.0370 | 0.0049 | 0.0081 | 0.0140 | 0.0137
ASOT-$k$ | 0.0455 | 0.0491 | 0.0863 | 0.0056 | 0.0256 | 0.0288 | 0.0053
eASOT-$k$ | 0.0343 | 0.0131 | 0.0587 | 0.0018 | 0.0075 | 0.0215 | 0.0039
eASOT-ML | 0.0087 | 0.0271 | 0.0316 | 0.0011 | 0.0017 | 0.0004 | 0.0052
eASOT-DL | 0.0045 | 0.0161 | 0.0146 | 0.0074 | 0.0087 | 0.0178 | 0.0040

Table 6: Average training time cost of graph datasets (sec).

Methods | MUTAG | PTC-MR | ENZYMES | BZR | COX2 | NCI1 | ZINC
---|---|---|---|---|---|---|---
ASOT-ML | 1.98 | 2.59 | 1.52 | 8.01 | 12.83 | 3.54 | 1.95
ASOT-DL | 5.97 | 6.07 | 5.86 | 5.12 | 4.78 | 7.09 | 0.41
ASOT-$k$ | 0.08 | 0.10 | 0.18 | 0.20 | 0.21 | 0.87 | 33.22

## References

  * [1] Cédric Villani.  Optimal transport: Old and new.  Springer, New York, 2008.
  * [2] Gabriel Peyré and Marco Cuturi.  Computational optimal transport.  Foundations and Trends in Machine Learning, 11(5-6):355–607, 2019\.
  * [3] Martin Arjovsky, Soumith Chintala, and Léon Bottou.  Wasserstein generative adversarial networks.  In International Conference on Machine Learning (ICML), 2017.
  * [4] Aude Genevay, Gabriel Peyré, and Marco Cuturi.  Learning generative models with Sinkhorn divergences.  In International Conference on Artificial Intelligence and Statistics (AISTATS), 2018.
  * [5] Matteo Togninalli, Elisabetta Ghisu, Felipe Llinares-López, Bastian Rieck, and Karsten Borgwardt.  Wasserstein Weisfeiler-Lehman graph kernels.  In Annual Conference on Neural Information Processing Systems (NeurIPS), 2019.
  * [6] Soheil Kolouri, Navid Naderializadeh, Gustavo K. Rohde, and Heiko Hoffmann.  Wasserstein embedding for graph learning.  In International Conference on Learning Representations (ICLR), 2021\.
  * [7] Bing Su and Gang Hua.  Order-preserving Wasserstein distance for sequence matching.  In Proceedings of the IEEE conference on computer vision and pattern recognition (CVPR), 2017.
  * [8] Mitsuhiko Horie and Hiroyuki Kasai.  Auto-weighted sequential Wasserstein distance and application to sequence matching.  In European Signal Processing Conference (EUSIPCO), 2022.
  * [9] 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.
  * [10] Soheil Kolouri, Kimia Nadjahi, Umut Simsekli, Roland Badeau, and Gustavo Rohde.  Generalized sliced Wasserstein distances.  In Annual Conference on Neural Information Processing Systems (NeurIPS), 2019.
  * [11] Zhongxi Fang, Sun Xun Huang, Jianming, and Hiroyuki Kasai.  Wasserstein graph distance based on L1-approximated tree edit distance between Weisfeiler-Lehman subtrees.  In AAAI Conference on Artificial Intelligence (AAAI), 2023.
  * [12] Jianming Huang, Zhongxi Fang, and Hiroyuki Kasai.  LCS graph kernel based on Wasserstein distance in longest common subsequence metric space.  Digital Signal Processing, 189:108281, 2021.
  * [13] Jianbo Ye, Panruo Wu, James Z Wang, and Jia Li.  Fast discrete distribution clustering using Wasserstein barycenter with sparse support.  IEEE Transactions on Signal Processing, 65(9):2317–2332, 2017.
  * [14] Takumi Fukunaga and Hiroyuki Kasai.  Wasserstein $k$-means with sparse simplex projection.  In International Conference on Pattern Recognition (ICPR), 2020\.
  * [15] Chi-Heng Lin, Mehdi Azabou, and Eva L Dyer.  Making transport more robust and interpretable by moving data through a small number of anchor points.  Proceedings of machine learning research, 139:6631, 2021.
  * [16] François-Pierre Paty and Marco Cuturi.  Subspace robust Wasserstein distances.  In International Conference on Machine Learning (ICML), 2019.
  * [17] Richard Sinkhorn and Paul Knopp.  Concerning nonnegative matrices and doubly stochastic matrices.  Pacific Journal of Mathematics, 21(2):343–348, 1967.
  * [18] Marco Cuturi.  Sinkhorn distances: Lightspeed computation of optimal transport.  In Annual Conference on Neural Information Processing Systems (NeurIPS), 2013.
  * [19] 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):1111–A1138, 2015.
  * [20] Jason Altschuler, Jonathan Niles-Weed, and Philippe Rigollet.  Near-linear time approximation algorithms for optimal transport via Sinkhorn iteration.  In Annual Conference on Neural Information Processing Systems (NeurIPS), 2017.
  * [21] Matthias Fey and Jan Eric Lenssen.  Fast graph representation learning with pytorch geometric.  arXiv preprint arXiv:1903.02428, 2019.
  * [22] The Popular Educator, volume 4, page 196.  John Cassel, 1855.
  * [23] Meyer Scetbon, Michael Elad, and Peyman Milanfar.  Deep K-SVD denoising.  IEEE Transactions on Image Processing, 30:5944–5955, 2021.
  * [24] Ingrid Daubechies, Michel Defrise, and Christine De Mol.  An iterative thresholding algorithm for linear inverse problems with a sparsity constraint.  Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 57(11):1413–1457, 2004.
  * [25] Asim Kumar Debnath, Rosa L Lopez de Compadre, Gargi Debnath, Alan J Shusterman, and Corwin Hansch.  Structure-activity relationship of mutagenic aromatic and heteroaromatic nitro compounds. correlation with molecular orbital energies and hydrophobicity.  Journal of medicinal chemistry, 34(2):786–797, 1991.
  * [26] Christoph Helma, Ross D. King, Stefan Kramer, and Ashwin Srinivasan.  The predictive toxicology challenge 2000–2001.  Bioinformatics, 17(1):107–108, 2001.
  * [27] Ida Schomburg, Antje Chang, Christian Ebeling, Marion Gremse, Christian Heldt, Gregor Huhn, and Dietmar Schomburg.  BRENDA, the enzyme database: updates and major new developments.  Nucleic acids research, 32(suppl_1):D431–D433, 2004.
  * [28] Jeffrey J Sutherland, Lee A O’brien, and Donald F Weaver.  Spline-fitting with a genetic algorithm: A method for developing classification structure-activity relationships.  Journal of chemical information and computer sciences, 43(6):1906–1915, 2003.
  * [29] Nikil Wale, Ian A Watson, and George Karypis.  Comparison of descriptor spaces for chemical compound retrieval and classification.  Knowledge and Information Systems, 14:347–375, 2008.
  * [30] Xavier Bresson and Thomas Laurent.  A two-step graph convolutional decoder for molecule generation.  arXiv preprint arXiv:1906.03412, 2019.
  * [31] Christopher Morris, Nils M. Kriege, Franka Bause, Kristian Kersting, Petra Mutzel, and Marion Neumann.  Tudataset: A collection of benchmark datasets for learning with graphs.  In International Conference on Machine Learning (ICML), 2020.
  * [32] Yann LeCun.  The MNIST database of handwritten digits.  http://yann. lecun. com/exdb/mnist/, 1998.
  * [33] Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka.  How powerful are graph neural networks?  In International Conference on Learning Representations (ICLR), 2019\.
  * [34] Nicolas Bonneel, Michiel Van De Panne, Sylvain Paris, and Wolfgang Heidrich.  Displacement interpolation using lagrangian mass transport.  In Proceedings of the Special Interest Group on Computer Graphics and Interactive Techniques (SIGGRAPH), 2011.
  * [35] Neal Parikh, Stephen Boyd, et al.  Proximal algorithms.  Foundations and Trends in Optimization, 1(3):127–239, 2014.
  * [36] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In International Conference on Machine Learning (ICML), 2015.
  * [37] David Sculley.  Web-scale k-means clustering.  In International Conference on World wide web (WWW), 2010.