Recovery Bounds on Class-Based Optimal Transport: A Sum-of-Norms Regularization Framework

TL;DR

This paper introduces a new convex optimal transport framework with sum-of-norms regularization that effectively recovers class structures, offers scalable algorithms, and enhances robustness and data structure preservation.

Contribution

It proposes a novel sum-of-norms regularized convex OT model with theoretical recovery guarantees and an efficient proximal algorithm.

Findings

Improved class structure preservation in data

Enhanced robustness to data geometry

Scalable computation of OT plans

Abstract

We develop a novel theoretical framework for understating OT schemes respecting a class structure. For this purpose, we propose a convex OT program with a sum-of-norms regularization term, which provably recovers the underlying class structure under geometric assumptions. Furthermore, we derive an accelerated proximal algorithm with a closed-form projection and proximal operator scheme, thereby affording a more scalable algorithm for computing optimal transport plans. We provide a novel argument for the uniqueness of the optimum even in the absence of strong convexity. Our experiments show that the new regularizer not only results in a better preservation of the class structure in the data but also yields additional robustness to the data geometry, compared to previous regularizers.

Tables

Table 1. Table 1: Accuracy scores on MNIST and USPS.

model	real1	real2	real3	real4
OT-SON	0.550	0.564	0.608	0.628
OT-l1l2	0.421	0.507	0.500	0.621
OT-lpl1	0.457	0.521	0.516	0.592
OT-Sinkhorn	0.414	0.521	0.508	0.621

Table 2. Table 2: Accuracy scores on Caltech Office (S: source, T: target).

S$\to$T	OT-SON	OT-l1l2	OT-lpl1	OT-Sinkhorn
A$\to$C	0.3552	0.3229	0.3565	0.3018
A$\to$D	0.3274	0.3097	0.3539	0.3008
A$\to$W	0.3144	0.2577	0.3092	0.2422
C$\to$A	0.4601	0.3714	0.4120	0.3548
C$\to$D	0.3982	0.3274	0.4424	0.3362
C$\to$W	0.3556	0.2474	0.3144	0.2474
D$\to$A	0.3248	0.2812	0.3248	0.2812
D$\to$C	0.2720	0.2658	0.2956	0.2621
D$\to$W	0.7216	0.7628	0.5876	0.7525
W$\to$A	0.2661	0.2225	0.2616	0.2210
W$\to$C	0.2149	0.1962	0.2124	0.2012
W$\to$D	0.8230	0.7964	0.6637	0.8141

Equations

$$\mu =\sum _{\alpha =1}^K{p}_{\alpha }{\mu }_{\alpha },\nu =\sum _{\beta =1}^K{q}_{\beta }{\nu }_{\beta }$$

$$\underset{\mathbf{X}\in B(\mathbf{\mu },\mathbf{\nu })}{\min }⟨\mathbf{D},\mathbf{X}⟩.$$

$${\mathbf{X}}^{\ast }=\arg \underset{\mathbf{X}\in B(\mathbf{\mu },\mathbf{\nu })}{\min }⟨\mathbf{D},\mathbf{X}⟩$$

$$+\lambda \left(\sum _{l,k}{R}_{l,k}\Vert {\mathbf{x}}_l-{\mathbf{x}}_k{\Vert }_2+\sum _{l,k}{S}_{l,k}\Vert {\mathbf{x}}^l-{\mathbf{x}}^k{\Vert }_2\right),$$

$$\frac{D-d-2\omega }{K^{\frac{3}{2}}}\ge C\omega ,$$

$$\frac{D^2-d^2}{K^{\frac{3}{2}}}\ge C\omega \sqrt{(E+\omega {)}^2+1}\log (nK),$$

$$\sum _{l=1}^k{D}_{{\alpha }_l{\alpha }_{l+1}}>\sum _{l=1}^k{D}_{{\alpha }_l{\alpha }_l}+k\delta .$$

$${\Lambda }_{\alpha ,\beta }={\left(T_{\alpha ,\beta }+\frac{{\omega }_{\alpha }+R{\omega }_{\beta }}{\sqrt{{\omega }_{\beta }^2+{\omega }_{\beta }^2}}\right)}^{-1},$$

$${\varphi }_{\rho ,\mathbf{\zeta },\mathbf{\eta }}(\mathbf{p},\mathbf{q})=⟨\mathbf{p},\mathbf{\zeta }⟩+⟨\mathbf{q},\mathbf{\eta }⟩+\rho \Vert \mathbf{p}-\mathbf{q}{\Vert }_2,$$

$$\sum _{l\ne k}{\varphi }_{R_{l,k},\frac{1}{2(n-1)}{\mathbf{d}}_l,\frac{1}{2(n-1)}{\mathbf{d}}_k}({\mathbf{x}}_l,{\mathbf{x}}_k)+$$

$$\sum _{l\ne k}{\varphi }_{S_{l,k},\frac{1}{2(m-1)}{\mathbf{d}}^l,\frac{1}{2(m-1)}{\mathbf{d}}^k}({\mathbf{x}}^l,{\mathbf{x}}^k),$$

$$S^{(n)}=\left\{\mathbf{x}=(x_i\ge 0{)}_{i=1}^n\mid \sum _i{x}_i=1\right\},$$

$$\underset{x\in {\mathbb{R}}^D}{\min }\sum _{p=1}^P{\varphi }_p(x)\text{s}tx\in \bigcap _{q=1}^Q{S}_q,$$

$$\mathbf{x}_{t+1}=\left\{\begin{array}[]{lc}\mathrm{prox}_{\mu\phi_{p_{t}}}\left(\mathbf{x}_{t}+\mu\mathbf{g}_{p_{t}}\right),&\phi_{p_{t}}\ \text{is selected}\\ \mathrm{proj}_{S_{q_{t}}}\left(\mathbf{x}_{t}+\mu\mathbf{h}_{q_{t}}\right)&S_{q_{t}}\ \text{is selected}\end{array}\right.,$$

$${\mathbf{a}}_t=\rho \frac{{\mathbf{x}}_t-{\mathbf{x}}_{t+1}}{\mu }-\alpha \left(\sum _n{\mathbf{g}}_n+\sum _m{\mathbf{h}}_m\right),$$

$$0\in \sum _{p=1}^P\partial {\varphi }_p(\mathbf{x})+\sum _{q=1}^Q\partial I_{S_q}(\mathbf{x})$$

$$\sum _p{\Vert {\mathbf{g}}_p^{\ast }\Vert }^2+\sum _q{\Vert {\mathbf{h}}_q^{\ast }\Vert }^2=O({\beta }^{2}R)$$

$$\mathbb{E}\left[\sum _p{\varphi }_p({\bar{\mathbf{x}}}_t)\right]-{\varphi }^{\ast }\le c\beta P\left(\frac{\eta }{t}+\sqrt{\frac{\beta \mu \eta }{t}}\right),$$

$$\mathbb{E}\left[\sum _q{\mathrm{dist}}^2({\bar{\mathbf{x}}}_t,S_q)\right]\le c\beta \mu P\frac{\eta }{t},$$

$${\mathcal{T}}_{\lambda }(\mathbf{a},\mathbf{b})=\left(\frac{\mathbf{a}+\mathbf{b}}{2}+{\mathcal{T}}_{\lambda }\left(\frac{\mathbf{a}-\mathbf{b}}{2}\right),\frac{\mathbf{a}+\mathbf{b}}{2}-{\mathcal{T}}_{\lambda }\left(\frac{\mathbf{a}-\mathbf{b}}{2}\right)\right),$$

$$\underset{\mathbf{X}\in {\mathbb{R}}_{\ge 0}^{n\times n}}{\min }⟨\mathbf{D},\mathbf{X}⟩+$$

$$\lambda \left(\sum _{i,i^{\prime }}{R}_{i,i^{\prime }}\Vert {\mathbf{x}}_i-{\mathbf{x}}_{i^{\prime }}{\Vert }_2+\sum _{j,j^{\prime }}{S}_{j,j^{\prime }}\Vert {\mathbf{x}}^j-{\mathbf{x}}^{j^{\prime }}{\Vert }_2\right)$$

$$+\frac{\theta }{2}\left({\Vert X1-\mathbf{\mu }\Vert }_2^2+{\Vert {\mathbf{X}}^{T}1-\mathbf{\nu }\Vert }_2^2\right)$$

$$\underset{X_{\alpha ,\beta }\ge 0}{\min }\sum _{\alpha ,\beta }{n}_{\alpha }{m}_{\beta }{X}_{\alpha ,\beta }{D}_{\alpha ,\beta }+$$

$$\lambda \left(\sum _{\alpha ,{\alpha }^{\prime }}{R}_{\alpha ,{\alpha }^{\prime }}\Vert {\mathbf{x}}_{\alpha }-{\mathbf{x}}_{{\alpha }^{\prime }}{\Vert }_M+\sum _{\beta ,{\beta }^{\prime }}{S}_{\beta ,{\beta }^{\prime }}\Vert {\mathbf{x}}^{\beta }-{\mathbf{x}}^{{\beta }^{\prime }}{\Vert }_N\right)$$

$$+\frac{\theta }{2}\left(\sum _{\alpha }{n}_{\alpha }{\left({\mathbf{a}}_M^T{\mathbf{x}}_{\alpha }-{\mu }_{\alpha }\right)}^2+\sum _{\beta }{m}_{\beta }{\left({\mathbf{a}}_N^T{\mathbf{x}}^{\beta }-{\nu }^{\beta }\right)}^2\right)$$

$$R_{\alpha ,{\alpha }^{\prime }}=\sum _{i\in C_{\alpha },i^{\prime }\in C_{{\alpha }^{\prime }}}{R}_{i,i^{\prime }},S_{\beta ,{\beta }^{\prime }}=\sum _{j\in D_{\beta },j^{\prime }\in D_{{\beta }^{\prime }}}{S}_{j,j^{\prime }}$$

$$D_{\alpha ,\beta }=\frac{\sum _{i\in C_{\alpha },j\in D_{\beta }}{D}_{i,j}}{n_{\alpha }{m}_{\beta }},{\mu }_{\alpha }=\frac{\sum _{i\in C_{\alpha }}{\mu }_i}{n_{\alpha }},{\nu }^{\beta }=\frac{\sum _{j\in D_{\beta }}{\nu }_j}{m_{\beta }},$$

$$\underset{Y_{\alpha ,\beta }\ge 0}{\min }\sum _{\alpha ,\beta }{Y}_{\alpha ,\beta }{D}_{\alpha ,\beta }$$

$$s.t$$

$$q_{\beta }:1^T{\mathbf{y}}^{\beta }={\sigma }^{\beta },p_{\alpha }:1^T{\mathbf{y}}_{\alpha }={\sigma }_{\alpha }$$

$$R_{i,\alpha }=\sum _{i^{\prime }\in C_{\alpha }}{R}_{i,i^{\prime }},S_{j,\beta }=\sum _{j^{\prime }\in D_{\beta }}{S}_{j,j^{\prime }}$$

$$\underset{X_{\alpha ,\beta }\ge 0}{\min }\sum _{\alpha ,\beta }{n}_{\alpha }{m}_{\beta }{X}_{\alpha ,\beta }{D}_{\alpha ,\beta }+$$

$$\lambda \left(\sum _{\alpha ,{\alpha }^{\prime }}{R}_{\alpha ,{\alpha }^{\prime }}\Vert {\mathbf{x}}_{\alpha }-{\mathbf{x}}_{{\alpha }^{\prime }}{\Vert }_M+\sum _{\beta ,{\beta }^{\prime }}{S}_{\beta ,{\beta }^{\prime }}\Vert {\mathbf{x}}^{\beta }-{\mathbf{x}}^{{\beta }^{\prime }}{\Vert }_N\right)$$

$${\mathbf{a}}_M^T{\mathbf{x}}_{\alpha }={\mu }_{\alpha },{\mathbf{a}}_N^T{\mathbf{x}}^{\beta }={\nu }^{\beta }$$

$$T_{\alpha }=\frac{1}{n_{\alpha }{m}_{\pi (\alpha )}}\sum _{{\alpha }^{\prime }\ne \alpha }\frac{R_{\alpha ,{\alpha }^{\prime }}\frac{{\sigma }_{\alpha }}{n_{\alpha }{m}_{\pi (\alpha )}}}{\sqrt{{\left(\frac{{\sigma }_{\alpha }}{n_{\alpha }{m}_{\pi (\alpha )}}\right)}^2+{\left(\frac{{\sigma }_{{\alpha }^{\prime }}}{n_{{\alpha }^{\prime }}{m}_{\pi ({\alpha }^{\prime })}}\right)}^2}},$$

Full text

Machine Learning, ICML

1 Introduction

Optimal transport (OT) is a classical mathematical discipline for discovering a transport map from a source distribution to a target distribution with a minimum cost of the transport. It has recently been successfully used in various applications in computer vision, texture analysis, tomographic reconstruction and clustering, as documented in the recent surveys (Kolouri et al.,, 2017) and (Solomon,, 2018). In many of these applications, OT exploits the geometry of the underlying spaces to effectively yield improved performance over the alternative of obviating it.

The main purpose of this paper is to provide a theoretical foundation for OT in such geometrically-aware conditions. We focus on a scenario where a common, potentially hidden class structure is present in both domains. Examples are abundant, such as (Long et al.,, 2022) and (Ott et al.,, 2022), where respecting the class structure can significantly improve the performance. In our setup, classes are associated with well separated regions of the data space, called components, on which the source and target distributions are supported. Each component is associated with a class and hence there is a correspondence between the components of the source and target domains. Such a model may become relevant after a suitable re-representation (embedding) of the data in a latent space.

Our central question is to identify conditions on the geometry of the components, under which OT can be performed in a polynomial time, with an additional property that each source sample is mapped to its corresponding component in the target domain. For this purpose, we introduce and study a two-stage procedure, where in the first stage the components and their association is recovered and in the second one, OT is exclusively performed over the corresponding component pairs. Any off-the-shelf OT method can be used in the second stage, and hence our results mainly pertain to the recovery condition in the first stage. For this, we introduce a novel convex optimization framework, combining two popular procedures: the celebrated Kantorovich relaxation scheme of OT and the well-known sum-of-norms (SON) formulation for vector clustering (Lindsten et al., 2011a, ; Hocking et al.,, 2011), used as a regularization. Accordingly, we show for the first time that sufficiently well-separated and associated components in the two domains can be recovered by a regularized OT scheme, which naturally enjoys a polynomial algorithm due to convexity. No such results, to the best of our knowledge, are known for other regularizers. We also experimentally show that our regularizer does not only yield a better class structure preservation, but also provides additional robustness compared to other class-based regularizers

**Computational benefits: ** Despite a convex nature, the improvements of OT often come at a significant computational cost. We further argue that the SON regularization of OT also enjoys practical computational benefits by proposing a novel stochastic incremental algorithm. First, we construct an abstract stochastic framework that is based on a combination of proximal and projection iterations, for which we give a generic proof of convergence at rate $O(1/T)$. Subsequently, we specialize this general scheme for our problem, which leads to an algorithm with computationally low-cost iterations. Beyond the proposed SON-based framework, our proximal scheme can be used to avoid the reported convergence difficulties of gradient-based methods (Patrascu and Necoara,, 2018).

Summary of contributions: Our main contributions can be summarized as follows:

i. In section 2.1, we develop a theoretical framework for class-based OT, where we introduce the concept of multi-class recovery schemes.

ii. As an instance of a multi-class recovery scheme, we propose in Section 2.2 a new regularized formulation of OT that recovers a class structure typically arising in real-world problems.

iii. We derive the first rigorous results for recovering an OT plan that respects class structure, presented in section 3, with more details provided in the appendix (section A).

iv. We develop in Section 4 a general accelerated stochastic incremental proximal-projection optimization scheme, for which we give a proof of convergence at a rate $O(1/T)$ without a decaying step size. We specialize the general scheme with an explicit closed form of proximal operators and fast projections to yield a scalable stochastic incremental algorithm for computing our OT formulation.

v. In section 5 and further in the appendix (section C), we investigate the algorithm on several synthetic and benchmark data sets, and demonstrate the benefits of the new regularizer.

vi. In the appendix (section D), we develop a new proof for the uniqueness of the optimal solution of our convex formulation in spite of its non-strong convexity, with wide applicability in other model recovery studies.

1.1 Relation to Literature

OT, first proposed by Monge as an analysis problem (Monge,, 1781), has become a classic topic in probability and statistics. A comprehensive introduction and theoretical framework can be found in (Villani,, 2008; Santambrogio,, 2015; Peyré et al.,, 2019). Theoretical works on extensions of OT have only recently received more attention. (Redko et al.,, 2017), for example, provides an analysis of OT in the context of domain adaptation. In (Gordaliza et al.,, 2019), OT is studied from the fairness point of view. To the best of our knowledge, there is no study on the geometrically-aware extensions (regularizations) of OT, as this paper offers.

Regularized OT: A case for exploring new regularizers was made in (Courty et al.,, 2017) in the context of domain adaptation applications. In (Blondel et al.,, 2018) the primal and dual formulations of OT are regularized with a strongly convex term, and the constraints are relaxed with smooth approximations. (Dessein et al.,, 2018) also propose a framework to solve discrete optimal transport problems with smooth convex regularization. Regularization is often introduced to promote sparsity. The L2 regularization for example yields a sparse plan. It has notably been used in a doubly stochastic scheme in (Seguy et al.,, 2018). However, these techniques are not tailored to the underlying class structure of the data, which may not be known in advance. The issue of a hidden class structure is addressed in (Courty et al.,, 2017) by the so-called Laplacian regularizer. However, no theoretical study is provided. We further illustrate the differences between the Laplacian regularizer and our framework in the experiments. It is also notable that our formulation differs from the partial domain adaption (Cao et al.,, 2018) and robust domain adaption (Balaji et al.,, 2020) settings developed to deal with outlier classes or samples. They do not generally take the class structure into account.

Computational aspects: Much attention has focused on efficient computational and numerical algorithms for OT, and a monograph focusing on this topic has appeared in (Peyré et al.,, 2019). In (Guo et al., 2020a, ), an ”accelerated primal-dual randomized coordinate descent (APDRCD)” algorithm is developed to solve the OT problem. An upper bound is also provided for the complexity of the algorithm and it is shown that it could be used for large-scale purposes. In (Courty et al.,, 2017), a generalized conditional gradient method is used to compute OT with the help of a couple of regularizers. Most notably Cuturi introduced an entropic regularizer and showed that its adoption with the Sinkhorn algorithm yields a fast computation of OT (Cuturi,, 2013); a theoretical guarantee that the Sinkhorn iteration computes the approximation in near linear time was also provided by (Altschuler et al.,, 2017). Screenkhorn algorithm proposed in (Alaya et al.,, 2019) performs screening to eliminate (and accelerate) the solution of the Sinkhorn algorithm. Another computational breakthrough was achieved by (Genevay et al.,, 2016) who gave a stochastic incremental algorithm to solve the entropic regularized OT problem. Proximal method has been used in OT but not in a stochastic scheme. Examples include (Alvarez-Melis et al.,, 2018) that uses notions of sub-modularity and also (Papadakis et al.,, 2014). Unlike the above works, our algorithm is tailored to the class-based OT framework, and exploits proximal/ projection operators in a stochastic way.

2 Problem Formulation and SON-regularized Optimal Transport

2.1 Problem Formulation

Consider two probability measures $\mu,\nu$ defined on data spaces ${\mathcal{Y}}^{s},{\mathcal{Y}}^{t}$, respectively referred to as the source and the target domains. Further, assume a non-negative function $d:{\mathcal{Y}}^{s}\times{\mathcal{Y}}^{t}\to\mathbb{R}_{\geq 0}$ evaluating the transport cost between the two domains. Classical Monge problem seeks a measurable transport map $T:{\mathcal{Y}}^{s}\to{\mathcal{Y}}^{t}$ such that the push-forward measure $T\sharp\mu$ of $\mu$ under $T$ coincides with $\nu$ and the expected transport cost $\mathbb{E}\left[d(Y,T(Y))\right]$ is minimized, where $Y\sim\mu$.

We similarly define the class-based OT problem:

Definition 1**.**

A $K$-class structure is a pair of mixtures of $K$ probability measures:

[TABLE]

on ${\mathcal{Y}}^{s}$ and ${\mathcal{Y}}^{t}$, respectively, with a one-to-one correspondence $\pi$ between the components of the two domains (i.e. $\mu_{\alpha},\nu_{\beta}$). We define a solution to the $K-$class Monge problem as a transport map $T:{\mathcal{Y}}^{s}\to{\mathcal{Y}}^{t}$ such that:

For every corresponding pair $(\mu_{\alpha},\nu_{\beta})\in\pi$ we have $T\sharp\mu_{\alpha}=\nu_{\beta}$, i.e $T$ transports any component to its corresponding component. 2. 2.

The expected transport cost $\mathbb{E}\left[d(Y,T(Y))\right]$ is minimized, where $Y\sim\mu$.

Note that a solution $T$ of the multi-class Monge problem depends on the components and their association. We consider a case that these are not provided:

Definition 2**.**

Given a transport cost function $d$, a Monge recovery scheme assigns to each pair $(\mu,\nu)$ of source and target distributions, a $K-$class structure $\left((\mu_{\alpha},p_{\alpha})_{\alpha},(\nu_{\beta},q_{\beta})_{\beta},\pi\right)$ for some $K\geq 0$, and its corresponding $K-$class Monge solution $T$ such that equation 2.1 holds true. A recovery scheme is said to recover a family ${\mathcal{F}}$ of multi-class structures if it assigns to any $(\mu,\nu)$ of ${\mathcal{F}}$ its corresponding mixture components and their association in ${\mathcal{F}}$. This naturally requires the components and their association to be unique in ${\mathcal{F}}$.

The main purpose of this paper is to introduce a wide family ${\mathcal{F}}$ of mixture models and a particular Monge recovery scheme recovering it in a polynomial time with the problem size.

2.1.1 Disjoint Supports and Two-Stage Schemes

Our focus is on mixture models with disjoint supports. For this case, we make the following straightforward observation:

Proposition 1*.*

Consider a multi-class structure $\left((\mu_{\alpha},p_{\alpha})_{\alpha},(\nu_{\beta},q_{\beta})_{\beta},\pi\right)$ and suppose that the supports $\mathscr{S}_{\alpha}$ of $\mu_{\alpha}$ are disjoint, and the same property holds for the supports $\mathscr{T}_{\beta}$ of $\nu_{\beta}$. Then the multi-class Monge solution comprises the restricted maps $T_{\alpha}:{\mathcal{S}}_{\alpha}\to{\mathcal{T}}_{\beta}$, for $(\alpha,\beta)\in\pi$, each being the solution of the conventional Monge problem for $(\mu_{\alpha},\nu_{\beta})$.

For such families of structures, the above observation suggests a two-stage recovery scheme: Given a pair of distributions $(\mu,\nu)$, first recover the components $(\mu_{\alpha}),(\nu_{\beta})$ and their association $\pi$. Next, solve the conventional Monge problem over each associated pair. We adopt this strategy, and as a rich theory already exists for the conventional Monge problem, we mainly focus on the first stage, i.e. recovering the components and their association.

2.2 Optimal Transport with SON

A major challenge in optimal transport is that only a finite number of samples from each distribution is given. Otherwise the distributions are unknown. Consider two finite sets $\{\mathbf{y}_{i}^{s}\}_{i=1}^{m},\{\mathbf{y}_{j}^{t}\}_{j=1}^{n}$ of points, respectively sampled from $\mu$ and $\nu$. Let $\mathbf{D}=(D_{ij}=d(\mathbf{y}^{s}_{i},\mathbf{y}^{t}_{j}))$ be the $m\times n$ matrix of transport costs. We denote the $i^{\textrm{th}}$ row and $j^{\textrm{th}}$ column of $\mathbf{D}$ by $\mathbf{d}^{i}$ and $\mathbf{d}_{j}$, respectively. We let the positive probability masses $\mu_{i},\nu_{j}$ be respectively assigned to the data points $\mathbf{y}_{i}^{s}$ and $\mathbf{y}^{t}_{j}$. In this discrete setup, the Monge problem is often solved with a linear programming relaxation scheme known as the Kantorovich problem:

[TABLE]

Here, the variable matrix $\mathbf{X}=(x_{i,j})$ is called the transport map and $B(\bm{\mu},\bm{\nu})=\{\mathbf{X}\in R^{m\times n},\mathbf{X}\mathbf{1}_{n^{s}}=\bm{\mu},\mathbf{X}^{T}\mathbf{1}_{n^{t}}=\bm{\nu}\}$ is the set of all discrete coupling distributions between $\bm{\mu}$ and $\bm{\nu}$, respectively denoting the vectors of elements $\mu_{i},\nu_{j}$. Moreover, $\langle\mathbf{D},\mathbf{X}\rangle=\mathrm{Tr}(\mathbf{D}^{T}\mathbf{X})=\sum_{i,j}X_{ij}D_{ij}$ is the Euclidean inner product of two matrices. In an ideal case, one hopes that the optimal solution for $\mathbf{X}$ become an assignment (permutation matrix) in which case it is seen to coincide with the solution of the Monge problem for the empirical distributions.

2.2.1 Multi-Class Recovery Scheme

Now, we introduce our multi-class recovery scheme by the following convex optimization problem:

[TABLE]

where $\mathbf{x}_{l}$ and $\mathbf{x}^{k}$ denote the (transpose of the) $l^{\textrm{th}}$ row and $k^{\textrm{th}}$ column of $\mathbf{X}$, respectively. Compared to equation 2.2, a regularization term with a tuning parameter $\lambda>0$ is introduced, known as sum of norms (SON). SON is well-known for its clustering properties and hence equation 2.2.1 combines the class discovery properties of SON with OT. The positive kernel coefficients $S_{l,k},R_{l,k}$ are also introduced to incorporate class prior information.

The effect of the SON regularization in equation 2.2.1 is explained in (Lindsten et al., 2011b, ; Panahi et al.,, 2017). In short, it enforces many vanishing regularization terms (sparsity), hence yielding identical columns and identical rows in the solution. In other words, the resulting map $\mathbf{X}^{*}$ after a suitable permutation of rows and columns is a block matrix with constant values in each block. The block structure reflects the discovered components in the source and target domains. Under suitable conditions , the constraints of the Kanturovich problem will further force many blocks to be zero. If each row and column will contain exactly one non-zero block (i.e. a block diagonal matrix under a correct order), the solution further reflects an assignment $\pi$ between the components. In this way, equation 2.2.1 performs the first stage in the two-stage multi-class recovery scheme. This is made precise and proved in section 3.

The regularization parameter $\lambda$ sets a desired balance between the cluster structure and the underlying transport problem. When $\lambda=0$, equation 2.2.1 reduces to equation 2.2. For large values of $\lambda$, the SON regularization dominates the result. In a typical situation, $\lambda\to\infty$ results in all data in each domain assigned to the same cluster, hence a trivial transformation between single clusters in each domain. Smaller values of $\lambda$ lead to a larger number of identified classes in the solution.

As mentioned, the kernel coefficients $R_{l,k},S_{l,k}$ are related to the prior knowledge of the components. For example if a perfect knowledge of the components in the source domain is available, we may set $R_{l,k}=0$ if $\mathbf{y}^{s}_{l}$ and $\mathbf{y}^{s}_{k}$ belong to different components (classes), otherwise set $R_{l,k}=k_{s}(\mathbf{y}^{s}_{l},\mathbf{y}^{s}_{k})$ for a suitable (differentiable) kernel $k_{s}$. On the target side where no class information is ordinarily provided, we may set $S_{l,k}=k_{t}(\mathbf{y}^{t}_{l},\mathbf{y}^{t}_{k})$ for a suitable kernel $k_{t}$ of choice.

3 Main Theoretical Results

Geometric threshold: To provide the main geometric intuitions of our analysis, we start by a simplified case with well-separated components. Next, we will present a more extensive study, which is also used for proving the simplified result:

Theorem 3.1**.**

Consider a mixture of $K$ components in each domain with $m$ samples from each component, hence a total of $n=Km$ samples in each domain. Suppose that $(\alpha,\alpha)\in\pi$ for all $\alpha$.

*Consider the $\ell_{2}$ distance, $d(\mathbf{y}_{1},\mathbf{y}_{2})=\|\mathbf{y}_{1}-\mathbf{y}_{2}\|_{2}$ and assume that the components are supported on spheres of radius $\omega$ and centered on $\bm{\theta}^{s}_{\alpha},\bm{\theta}^{t}_{\alpha}$ for $\alpha=1,2,\ldots,K$, in the source and target domains, respectively. With a suitable choice of $\lambda$, the solution of equation 2.2.1 recovers the *$K-$class structure if:

[TABLE]

for some universal constant $C$, where $D=\min_{\alpha\neq\beta}\|\bm{\theta}^{s}_{\alpha}-\bm{\theta}^{t}_{\beta}\|$ and $d=\max_{\alpha}\|\bm{\theta}^{s}_{\alpha}-\bm{\theta}^{t}_{\alpha}\|$. 2. 2.

*Consider the $\ell^{2}_{2}$ distance, $d(\mathbf{y}_{1},\mathbf{y}_{2})=\|\mathbf{y}_{1}-\mathbf{y}_{2}\|^{2}_{2}$. Assume that the components are isotropic Gaussian with means $\bm{\theta}^{s}_{\alpha},\bm{\theta}^{t}_{\alpha}$ for $\alpha=1,2,\ldots,K$, in the source and target domains, respectively. All variances are equal to $\omega^{2}$. With a probability higher than $1-\nicefrac{{1}}{{n^{10}}}$ the solution of equation 2.2.1 with a suitable choice of $\lambda$ recovers the *$K-$class structure if:

[TABLE]

where $E=\max_{\alpha,\beta}\|\bm{\theta}^{s}_{\alpha}-\bm{\theta}^{t}_{\beta}\|$.

Fig. 1 clarifies in a simple example the geometric meaning of the concepts used in the above result. As seen, the conditions in equation 3.2 and equation 3.1 require the associated clusters to be substantially closer to each other than the other clusters. As expected, there are differences between the $\ell_{2}$ and $\ell_{2}^{2}$ geometries. For example, the scale $E$ of the problem only appears in the latter. However, in both cases we may bound the maximum number of resolvable components $K$ by a purely geometric parameter that we refer to as the geometric threshold. In case 1), for example, defining $r=\frac{D-d-2\omega}{\omega}$, we observe that the geometric threshold is $O(r^{\frac{2}{3}})$. Part 2) has a similar interpretation.

Deterministic guarantee: Now, we present an extended deterministic result that is used to prove theorem 3.1. In many respects, the presented results are still simplified and the most comprehensive analysis is postponed to Appendix B. For simplicity, $K$ components with equal size (number of samples) $m$ is assumed. The resulting data partitions in the source and target domains are respectively denoted by $\{\mathscr{S}_{\alpha}\}$, $\{\mathscr{T}_{\beta}\}$. The total number of points in each domain is $n=mK$ . Further, $\mathscr{S}_{\alpha}$ is paired with $\{\mathscr{T}_{\alpha}\}$ for every $\alpha\in[K]$, i.e. $(\alpha,\alpha)\in\pi$. We investigate that the transport plan obtained by solving equation 2.2.1 consists of blocks, recovering both the sets of clusters $\{\mathscr{S}_{\alpha}\}$, $\{\mathscr{T}_{\beta}\}$ and their association. For this, we ensure that $X_{ij}$ remains zero for the $i^{\textrm{th}}$ data point in the source domain and $j^{\textrm{th}}$ data point in the target domain, belonging to unassociated clusters. Accordingly, we require the ideal solution to be the one with $X_{i,j}=X_{\alpha,\beta}$ for $i\in\mathscr{S}_{\alpha}$ and $j\in\mathscr{T}_{\beta}$, where $X_{\alpha,\beta}$ are constants satisfying $X_{\alpha,\beta}=0$ for $\beta\neq\alpha$.

For simplicity, we take $S_{j,j^{\prime}}=1$ everywhere and study two cases where $R_{i,i^{\prime}}=1$ holds true either everywhere (no kernel) or for $i,i^{\prime}$ belonging to the same cluster and $R_{i,i^{\prime}}=0$ otherwise (perfect kernels in the source domain). The general case is presented in Appendix B. Introducing an indicator variable $R$, the first case is referred to by $R=0$ and the second one by $R=1$. Note also that we assume the ideal solution to be feasible for the optimization problem in equation 2.2.1, which requires for every $i,i^{\prime}\in\mathscr{S}_{\alpha}$ and $j,j^{\prime}\in\mathscr{T}_{\alpha}$ that $\mu_{i}=\mu_{i^{\prime}}=\nu_{j}=\nu_{j^{\prime}}$. In Appendix B, we treat the general infeasible cases by considering a relaxation of equation 2.2.1.

In the context of recovery by the Kantorovich relaxation, a key concept is cyclical monotonicity (Villani,, 2008), which we slightly modify and state below:

Definition 3**.**

We say that a set of coefficients $D_{\alpha,\alpha^{\prime}}$ for $\alpha,\alpha^{\prime}\in[K]$ satisfies the $\delta-$strong cyclical monotonicity condition if for each simple loop $\alpha_{1}\to\alpha_{2}\to\ldots\to\alpha_{k}\to\alpha_{k+1}=\alpha_{1}$ with length $k>1$ we have

[TABLE]

Compared to the standard notion of cyclic monotonicity, we introduce a constant $\delta\geq 0$ on the right hand side of equation 3.3, which can be nonzero only when $(D_{\alpha,\beta})$ has a discrete or discontinuous nature. We apply this condition to the average distance of clusters given by $D_{\alpha,\beta}=\frac{1}{m^{2}}\sum\limits_{i\in\mathscr{S}_{\alpha},j\in\mathscr{T}_{\beta}}D_{i,j}.$

We denote by $\Delta$ the maximum of the values $\nicefrac{{\left\|\mathbf{d}_{i}-\mathbf{d}_{i^{\prime}}\right\|}}{{\sqrt{n}}}$ and $\nicefrac{{\left\|\mathbf{d}^{j}-\mathbf{d}^{j^{\prime}}\right\|}}{{\sqrt{n}}}$ where source points $i,i^{\prime}$ and target points $j,j^{\prime}$ belong to the same cluster and we remind that $\mathbf{d}_{i},\mathbf{d}^{j}$ respectively refer to the rows and columns of $\mathbf{D}$. We also define $\omega_{\alpha}:=\sum\limits_{i\in\mathscr{S}_{\alpha}}\mu_{i}=\sum\limits_{j\in D_{\alpha}}\nu_{j}$ and then take $T_{\alpha,\beta}=\sum\limits_{\gamma\in[K]}\left(\frac{R\omega_{\alpha}}{\sqrt{\omega_{\alpha}^{2}+\omega_{\gamma}^{2}}}+\frac{\omega_{\beta}}{\sqrt{\omega_{\beta}^{2}+\omega_{\gamma}^{2}}}\right)-\frac{1+R}{\sqrt{2}}$ . Finally, we define

[TABLE]

and take $\Lambda$ as its maximum over $\alpha\neq\beta$. Accordingly, we obtain the following result:

Theorem 3.2**.**

Suppose that $(D_{\alpha,\beta})$ is $\delta-$strongly cyclical monotone. Take $\lambda$ such that $\Delta\leq\lambda\sqrt{\nicefrac{{m}}{{K}}}$. Then, the solution of equation 2.2.1 is given by $X_{ij}=X_{\alpha,\beta}$ for $i\in\mathscr{S}_{\alpha}$ and $j\in\mathscr{T}_{\beta}$ satisfying one of the following two conditions:

We have $X_{\alpha,\beta}=\nicefrac{{\omega_{\alpha}}}{{m^{2}}}\delta_{\beta,\alpha}$ if $\Delta\sqrt{K}\leq\lambda\sqrt{m}\leq\Lambda\delta$ 2. 2.

Otherwise, we have $\delta\sum\limits_{\beta\neq\alpha}X_{\alpha,\beta}\leq\lambda(1+R)\sqrt{m}\sum\limits_{\alpha\neq\alpha^{\prime}}\sqrt{\omega^{2}_{\alpha}+\omega^{2}_{\alpha^{\prime}}}.$

Furthermore, the solution is unique in part 1 if all inequalities are strict.

Proof.

Proof can be found in section A.1. ∎

The first part of theorem 3.2 establishes ideal recovery. It can be understood in light of theorem 3.1, e.g., in part 1. The term $\Delta$ corresponds to $\omega$ in 3.1. It is the maximal cluster diameter in a geometry embedded by the vectors $\mathbf{d}_{i},\mathbf{d}^{j}$111Note that here $\mathbf{d}_{i},\mathbf{d}^{j}$ are treated as a special embedding of the source and target points, respectively, which may be called the distance embedding. As seen, we do not assume any inherent geometry in the two domains and instead rely on the induced $\ell_{2}$ geometry of the distance embedding.. On the other hand, $\delta$ corresponds to $D-d-2\omega$ which is the gain of the assignment. The second part gives an upper bound on the error $\sum\limits_{\beta\neq\alpha}X_{\alpha,\beta}$. Note that $\Delta$ is always smaller with $R=1$ compared to $R=0$, making the conditions less restrictive. This reflects the intuitive fact that introducing kernels simplifies the recovery.

4 Stochastic Incremental Algorithms

4.1 Accelerated Proximal-Projection Scheme

An important advantage of the framework in equation 2.2.1 is the possibility of applying stochastic optimization techniques. Since the objective term includes a large number of non-smooth SON terms, our stochastic optimization avoids calculating the (sub)gradient or the proximal operator of the entire objective function, which is numerically infeasible for large-scale problems. Our algorithm is obtained by introducing the following ”template function”:

[TABLE]

and noting that the objective function in equation 2.2.1 can be written as

[TABLE]

with a total number of $P=m(m-1)+n(n-1)$ summands in the form of the template function. This places the problem in the setting of finite sum optimization problems (Bottou et al.,, 2018). However, there are two obstacles to the application of stochastic optimization techniques: First, the terms in (LABEL:eq:fullobj) are not smooth, so gradient methods do not apply and second, equation 2.2.1 involves a fairly complex constraint. We address these issues in the following.

Non–smooth terms: We exploit the highly effective proximal methodology for optimizing non–smooth functions (Parikh and Boyd,, 2016; Combettes and Pesquet,, 2011) using a proximal operator. Defazio further gives a stochastic acceleration technique using proximal operators for unconstrained problems (Defazio,, 2016). In addition to its fast convergence, the main advantage of this scheme is its potential constant step size convergence in contrast to the ordinary stochastic gradient approach. It unfortunately does not address constrained optimization problems.

Constrained optimization: Facing a constrained optimization problem, the calculation of the proximal operators over the feasible set is numerically intractable. However, we observe an appealing structure in the constraint which lends itself to a more efficient stochastic implementation: Recalling the definition of an $n-$dimensional standard simplex

[TABLE]

we define the weighted cylinder-simplices $S_{l}(\mu)=\{\mathbf{X}\mid\mathbf{x}_{l}\in\mu S^{(n)}\}$ and $S^{k}(\nu)=\{\mathbf{X}\mid\mathbf{x}^{k}\in\nu S^{(m)}\}$ respectively corresponding to the $l^{\textrm{th}}$ row and $k^{\textrm{th}}$ column of $\mathbf{X}$ with weights $\mu,\nu\geq 0$. We then observe that the constraint set $B(\bm{\mu},\bm{\nu})$ is equal to $B(\bm{\mu},\bm{\nu})=(\bigcap_{l=1}^{m}S_{l}(\mu_{l}))\cap(\bigcap_{k=n}^{m}S^{k}(\nu_{k}))$, which is an intersection of $Q=m+n$ weighted cylinder-simplices.

In summary, the optimization problem in equation 2.2.1 can be written in the following abstract form:

[TABLE]

where each term $\phi_{p}$ denotes a template function term in the objective and each set $S_{q}$ is a weighted cylinder-simplex. The values of $P,Q$ in equation 2.2.1 are given above and $D=mn$. (Bertsekas,, 2011), (Wang and Bertsekas,, 2016) and (Patrascu and Necoara,, 2018) give general stochastic incremental schemes that combine gradient, proximal and projected schemes for optimizing such finite sum problems with convex constraints. These do not, however, use acceleration and their respective convergence is only guaranteed with a variable and vanishing step size, which is practically difficult to control and often yields extremely slow convergence.

Our proposed method: We herein jointly exploit the two ideas in (Defazio,, 2016) and (Wang and Bertsekas,, 2016) to obtain an accelerated proximal scheme for constrained framework in equation 2.2.1. Further, we shortly show in Lemma 4.2 that the proximal operator can be computed in closed form for our problem. Together with the projection to the simplex from (Condat,, 2016; Duchi et al.,, 2008), this gives a stochastic incremental algorithm with much less costly iterations.

We extend the acceleration techniques of unconstrained optimization as in the Defazio’s scheme (known as Point–SAGA) to the constrained setting. Point–SAGA utilizes individual ”memory” vectors for each term in the objective function, which store a calculated subgradient of a selected term in every iteration. These vectors are subsequently used as an estimate of the subgradient at subsequent iterations. We extend this scheme by introducing similar memory vectors to constraints. Each memory vector $\mathbf{h}_{m}$ for a constraint $S_{m}$ stores the last observed normal (separating) vector to $S_{m}$. At each iteration either an objective term $\phi_{p}$ or a constraint component $S_{q}$ is considered by random selection. Accordingly, we propose the following rule for updating the solution:

[TABLE]

where $t$ is the iteration number, $\mu>0$ is the fixed step size and $p_{t},q_{t}$ denote the selected index in this iteration (only one of them exists). At each iteration, the corresponding memory vector to the selected term is also updated. Depending on the choice of $\phi_{p_{t}}$ or $S_{q_{t}}$, either $\mathbf{g}_{p_{t}}\leftarrow\mathbf{g}_{p_{t}}+\mathbf{a}_{t}$ or $\mathbf{h}_{q_{t}}\leftarrow\mathbf{h}_{q_{t}}+\mathbf{a}_{t}$, where

[TABLE]

where $\rho\in(0\ 1)$ and $\alpha>0$ are design constants. The vector $\mathbf{a}_{t}$ consists of two parts: the first part $\rho\frac{\mathbf{x}_{t}-\mathbf{x}_{t+1}}{\mu}$ calculates a sub-gradient or a normal vector at point $\mathbf{x}_{t+1}$ corresponding to the selected term. The second term, the sum of the memory terms, implements acceleration. Our algorithm bears marked differences with Point-SAGA. While acceleration by the sum of memory vectors is also employed in Point-SAGA, it is moved in our scheme from the update rule of $\mathbf{x}_{t}$ to the update rule of $\mathbf{g}_{t}$. Also, the design parameters $\rho$ and $\alpha$ are introduced to improve convergence. Similar to Point-SAGA we only need to calculate the sum of memory terms once in the beginning and later update it by simple manipulations. As we later employ initialization of the memory vectors by zero, the first summation trivially leads to zero.

Convergence analysis:

We show that for a generic convex optimization problem of the form in equation 4.3, the algorithmic scheme in section 4.1 converges with a guaranteed rate, under the following mild assumptions:

Assumption 1*.*

The functions $\phi_{p}$ are $\beta-$Lipschitz.

Assumption 2*.*

We require the monotone inclusion problem

[TABLE]

to have a solution at $\mathbf{x}=\mathbf{x}^{*}$ with a finite optimal value $\phi^{*}$ and $\mathbf{g}_{p}^{*}\in\partial\phi_{p}(\mathbf{x}^{*})$ and $\mathbf{h}_{q}^{*}\in\partial I_{S_{q}}(\mathbf{x}^{*})$ satisfying $\sum_{p}\mathbf{g}_{p}^{*}+\sum_{q}\mathbf{h}_{q}^{*}=\mathbf{0}$. Furthermore, we assume that

[TABLE]

where $R=P+Q$ is the total number of terms.

Here, $\partial\phi(\mathbf{x})$ and $\partial I_{S}(\mathbf{x})$ respectively denote the subdifferential of the function $\phi$ and the cone of normal vectors to the set $S$ at $\mathbf{x}$. It is well-known that any solution to equation 4.6 is an optimal feasible solution to equation 4.3.

Assumption 3*.*

We assume that the algorithm is initialized with $\mathbf{g}_{p}=\mathbf{h}_{q}=\mathbf{0}$.

Then, we can show the following result.

Theorem 4.1**.**

Suppose that Assumption 1-3 are satisfied. Then for $\alpha,\rho,\mu>0$ and $\alpha<2(1-\rho)$ the following holds true:

Defining $\bar{\mathbf{x}}_{t}=\frac{1}{t}\sum\limits_{\tau=0}^{t-1}\mathbf{x}_{\tau}$ and $\eta=(1+\nicefrac{{Q}}{{P}})(\nicefrac{{\|\mathbf{x}_{0}-\mathbf{x}^{*}\|^{2}}}{{\beta\mu}}+\beta\mu R)$ we have

[TABLE]

where $c$ is a constant depending on $\rho,\alpha$ and the underlying constant in equation 4.7. 2. 2.

Moreover, $\sum\limits_{\tau=0}^{\infty}\mathbb{E}[\|\mathbf{x}_{\tau+1}-\mathbf{x}_{\tau}\|^{2}]\leq c\frac{\mu\beta P}{P+Q}\eta.$

Proof.

The proof is given in section B. ∎

Comment: As expected, the results only depend on $\beta\mu$, except for the optimality gap being linearly proportional to $\beta$. Applying this technique to our problem of interest and assuming that $m,n$ are of the same order, we observe that $P=O(n^{2})$ and $Q=O(n)$. Since many terms in our objective function are for regularization, it is fair to consider the relative optimality gap obtained by driving the optimality gap to the number of objective terms. We observe that this quantity is controlled by $\eta/t$. We conclude that $t\sim\eta$ iterations is required to achieve a desired relative optimality gap. The total feasibilty gap is controlled by $\beta\mu P\eta$. If we take $\beta\mu\sim\nicefrac{{1}}{{n}}$, we obtain $\eta\sim n$ and the relative optimality gap vanishes in $O(n)$ iterations. Then, the total optimality gap and feasibility gap will vanish in $O(n^{3})$ and $O(n^{2})$ iterations, respectively. In the absence of the regularization terms, we may reorganize the objective to have only $O(n)$ terms. In this case, taking $\beta\mu\sim\nicefrac{{1}}{{\sqrt{n}}}$, we get $\eta\sim O(\sqrt{n})$ and we require $O(n^{\frac{3}{2}})$ and $O(n)$ iterations to control the total optimality and feasibility gaps. Compared to the results of (Guo et al., 2020b, ), which establishes convergence in $O(n^{\frac{5}{2}})$ our convergence rates are better. Moreover, the $O(n^{2})$ dependence can be improved by reducing the number of terms in the objective function. It has been pointed out that in the sum of norms approach many terms may be redundant and only $O(n)$ terms corresponding to pre-selected pairs $(i,j)$ can be sufficient.

4.2 Proximal Operator for the SON-Regularized Kantorovich Relaxation

We next show that we can explicitly compute the proximal operator for each term in (LABEL:eq:fullobj):

Theorem 4.2**.**

The proximal operator of the template function $\phi_{\rho,\bm{\zeta},\bm{\eta}}$ is given by ${\mathcal{T}}_{\mu\rho}(\mathbf{p}-\mu\bm{\zeta},\mathbf{q}-\mu\bm{\eta}),$ where

[TABLE]

*and ${\mathcal{T}}_{\lambda}(\mathbf{c})$ is a thresholding operator given by $\frac{\|\mathbf{c}\|-\lambda}{\|\mathbf{c}\|}\mathbf{c}$ if $\|\mathbf{c}\|\geq\lambda$ and is zero otherwise. *

Proof.

Proof is found in section B.3. ∎

5 Experiments

We now investigate the effectiveness of our OT framework for the well-known domain adaptation problem which aims at improving the performance of a model in a target domain by utilizing the knowledge available in a different but related source domain. In the appendix (section C), we investigate the benefits of several other aspects of our framework, such as early stopping, class diversity, and unsupervised domain adaptation. We compare our method (OT-SON) with the other regularized optimal transport-based methods OT-l1l2, OT-lpl1 and OT-Sinkhorn, as developed and used in (Courty et al.,, 2017; Cuturi,, 2013; Perrot et al.,, 2016).

5.1 Domain Adaptation with Real-World Datasets

5.1.1 MNIST and USPS

In these experiments, we compare the different models on the real-world images of digits. For this, we consider the MNIST data as the source and the USPS data as the target. To further increase the difficulty of the problem, we use all $10$ classes of the source (MNIST) data, and we discard some of the classes of the target (USPS) data. In our experiments, each object (image) is represented by $256$ features. By discarding the different subsets from the USPS data, we consider several pairs of source and target datasets. i) real1: the USPS classes are $1,2,3,5,6,7,8$, ii) real2: the USPS classes are $0,2,4,5,6,7,9$, iii) real3: the USPS classes are $0,1,3,5,7,9$, and iv) real4: the USPS classes are: $0,1,3,4,6,8,9$. We note that these settings where the number of classes is different between the domains are the typical cases in practice. Therefore, our class-specific OT approach is more suitable and robust to class imbalance, as it avoids splitting a class in one domain among multiple classes in another domain.

The transformed source samples are used to train a 1-nearest neighbor classifier. We then use this (parameter-free) classifier to estimate the class labels in the target data and then compute the respective accuracy. Table 1 shows the best accuracy results for different OT-based models when using different values of the regularization parameters $\lambda_{1}$ and $\lambda_{2}$ (i.e., $\lambda_{1},\lambda_{2}\in\{10^{-5},...,10^{3}\}$)222For OT-l1l2 and OT-lpl1, $\lambda_{1}$ is the entropic regularization parameter and $\lambda_{2}$ is class regularization parameter.. Specifically, for each OT method we use different values of regularization parameters and we report the best accuracy achieved by that method. We observe, i) OT-SON yields the highest accuracy scores, and ii) it is significantly more robust to variation of the regularization parameters, in comparison to the other methods. Moreover, the other methods are prone to yielding numerical errors for small regularizations.

5.1.2 Caltech Office

A commonly used dataset for domain adaptation is an object recognition dataset known as Caltech Office (Saenko et al.,, 2010; Griffin et al.,, 2006) which consists of four different domains: A (Amazon, 958 samples) , W (Webcam, 295 samples), C (Caltech, 1123 samples), and D (DSLR, 157 samples). These domains have 10 classes of objects in common that are represented with two sets of features: SURF (800 features) and DeCAF (4096 features). Similar to the setting in (Courty et al.,, 2017), we use all possible pairs of the four domains as source and target with SURF features, and we remove the classes 3, 5, and 7 from the target domain. Similar to the previous study, we assume imbalanced source and target classes in order to make the task more realistic. Table 2 compares the classification accuracies achieved by our method with the scores obtained by the three other methods. We calculate the accuracy similar to Section 5.1.1. In these experiments we use class information in the source domain together with Gaussian kernels. Specifically, we set $R_{l,k}=0$ if $y_{l}^{s}$ and $y_{k}^{s}$ are not in the same classes and otherwise we set $R_{l,k}=\lambda\exp(-\|y_{l}^{s}-y_{k}^{s}\|^{2})$ for different values of $\lambda$. We observe that our method yields the best results in seven cases. In other cases, our method is still competitive compared to the alternatives. In addition, we conclude that the different methods may perform differently on different datasets. Even though the setting of imbalanced classes is more important in practice and we have focused more on that, for the sake of completeness, we also compare our method with the three other methods in a setting where the classes are balanced. We use the same number of classes in the source ($W$) and target ($C$) domains. The accuracy results for the different methods are respectively: i) OT-SON: 0.260, ii) OT-l1l2: 0.244, iii) OT-lpl1: 0.247, and iv) OT-Sinkhorn: 0.243. We observe that even in this setting, our method yields the highest accuracy.

6 Conclusion

We introduced a novel theoretical framework for OT in the presence of a multi-class structure in the two domains. Accordingly, we provided first theoretical guarantees for the recovery of the class structure, and developed constrained incremental algorithms which are generally suitable for non-smooth problems and enjoy theoretical convergence guarantees. Our experimental studies have substantiated the effectiveness of our proposed approach in different illustrative settings and datasets.

Acknowledgements

This work was partially supported by the Wallenberg AI, Autonomous Systems and Software Program (WASP) funded by the Knut and Alice Wallenberg Foundation. We would like to thank the anonymous reviewers for their constructive comments.

Appendix A Extension of Theorem 3.2

We consider the analysis of our proposed method for general kernel coefficient and cluster sizes. Hence, we respectively consider two partitions $\{C_{\alpha}\}$, $\{D_{\beta}\}$ of $[n],[m]$ with the same number of parts $K$. We denote the cardinalities of $C_{\alpha}$ and $D_{\beta}$ by $n_{\alpha}$ and $m_{\beta}$, respectively. Further, we consider a permutation $\pi$ on $[K]$ as the target of OT. Also, we address infeasibility by consider the following optimization:

[TABLE]

where $\theta>0$ is a design parameter and we remind that $\mathbf{x}_{i}=(X_{i,j})_{j}$, $\mathbf{x}^{j}=(X_{i,j})_{i}$, and $R_{i,i^{\prime}}$ and $S_{j,j^{\prime}}$ are positive kernel coefficients. Now, we introduce few intermediate optimizations to carry out the analysis. Define the following more general characteristic optimization:

[TABLE]

where

[TABLE]

[TABLE]

$\|\mathbf{x}\|_{O}=\sqrt{\mathbf{x}^{T}O\mathbf{x}}$, $N,M$ are diagonal matrices with $n_{\alpha}$ $m_{\beta}$ as diagonals, resepectively, $\mathbf{x}_{\alpha}=(X_{\alpha,\beta})_{\beta}$ and $\mathbf{x}^{\beta}=(X_{\alpha,\beta})_{\alpha}$, and $\mathbf{a}_{M}=(m_{\alpha})_{\alpha}$, $\mathbf{a}_{N}=(n_{\alpha})_{\alpha}$.

Further, define the ideal optimization:

[TABLE]

where $\sigma_{\alpha}=(n_{\alpha}\mu_{\alpha}+m_{\pi(\alpha)}\nu^{\pi(\alpha)})/2$, $\sigma^{\beta}=\sigma_{\pi^{-1}(\beta)}=(n_{\pi^{-1}(\beta)}\mu_{\pi^{-1}(\beta)}+m_{\beta}\nu^{\beta})/2$, and $\{p_{\alpha}\},\{q_{\beta}\}$ are dual variables. Also, define $\delta_{\alpha}=(\mu_{\alpha}n_{\alpha}-m_{\pi(\alpha)}\nu^{\pi(\alpha)})/2$, $\delta^{\beta}=-\delta_{\pi^{-1}(\beta)}=(m_{\beta}\nu^{\beta}-n_{\pi^{-1}(\beta)}\mu_{\pi^{-1}(\beta)})/2$ and $\bm{\delta}=(\delta_{\alpha})$. Finally, take

[TABLE]

In case $\bm{\delta}=0$, we may also define the tight characteristic optimization:

[TABLE]

which in this case, coincides with equation A when $\theta=\infty$. We also define

[TABLE]

Then, we have the following more general result:

Theorem A.1**.**

Suppose that $\tilde{D}_{\alpha,\alpha^{\prime}}=D_{\alpha,\pi(\alpha^{\prime})}$ satisfies the strong cyclical monotonicity condition, where for each simple loop $i_{1}\to i_{2}\to\ldots\to i_{k}\to i_{k+1}=i_{1}$ with length $k>1$ we have

[TABLE]

Then, the solution ${X_{\alpha,\beta}}$ of the characteristic optimization in equation A satisfies one of the following:

1. (a)

If $\bm{\delta}=0$ and $\theta=\infty$, we have $X_{\alpha,\beta}=\delta_{\beta,\pi(\alpha)}\frac{\sigma_{\alpha}}{m_{\alpha}n_{\beta}}$ with $\delta_{\ldotp,\ldotp}$ being the Kronecker index, if

[TABLE] 2. (b)

Otherwise,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

[TABLE] 2. 2.

The solution of equation A is given by $X_{ij}=X_{\alpha,\beta}$ if there exist positive constants $a,c,d$ such that $2a+c+d\leq 1$ and for all $i,i^{\prime}\in C_{\alpha}$ and $j,j^{\prime}\in D_{\beta}$,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Proof.

Denote the optimal value of equation A and equation A by $C_{0}$ and $C_{1}$, respectively. Also, notice that since $\tilde{D}_{\alpha,\alpha^{\prime}}$ satisfies the strong cyclical monotonicity condition, $Y_{\alpha,\beta}=\delta_{\beta,\pi(\alpha)}\sigma_{\alpha}$ is the solution of equation A and there exist dual variables $p_{\alpha},q_{\beta}$ such that

[TABLE]

Moreover,

[TABLE]

For part 1.a, we note that under the given conditions, the solution $X_{\alpha,\beta}$ of equation A coincides with that of equation A. Now, we show that $X^{\prime}_{\alpha,\beta}=\frac{Y_{\alpha,\beta}}{n_{\alpha}m_{\beta}}=\frac{\delta_{\beta,\pi(\alpha)}\sigma_{\alpha}}{n_{\alpha}m_{\beta}}$ satisfies with the dual parameters $p^{\prime}_{\alpha},q^{\prime}_{\beta}$, the optimality condition of equation A, which can be written as

[TABLE]

where $A_{\alpha,\beta}$ is the partial derivative at $X^{\prime}_{\alpha,\beta}$ of the SON term w.r.t $X_{\alpha,\beta}$. By direct calculation, we observe that

[TABLE]

It is now simple to check that under the given assumption, taking $p_{\alpha}^{\prime}=p_{\alpha}+\lambda T_{\alpha}$ and $q_{\beta}^{\prime}=q_{\beta}+\lambda U^{\beta}$ will satisfy the optimality conditions.

For part 1.b, we note that for the solution $X_{\alpha,\beta}$ of equation A,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $F(\ldotp)$ denotes the objective function in equation A. On the other hand for $X^{\prime}_{\alpha,\beta}=\frac{Y_{\alpha,\beta}}{n_{\alpha}m_{\beta}}=\frac{\delta_{\beta,\pi(\alpha)}\sigma_{\alpha}}{n_{\alpha}m_{\beta}}$, we have that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We conclude that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Lemma A.2 gives the result in part 1.

For part 2, notice that the optimality condition of $X_{\alpha,\beta}$ yields

[TABLE]

[TABLE]

where

[TABLE]

Define for $i,i^{\prime}\in C_{\alpha}$ and $j,j^{\prime}\in D_{\beta}$

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Also for $i\in C_{\alpha},i^{\prime}\in C_{\alpha^{\prime}}$ and $j\in D_{\beta},j^{\prime}\in D_{\beta^{\prime}}$, where $\alpha\neq\alpha^{\prime}$ and $\beta\neq\beta^{\prime}$, take $(\mathbf{z}_{ii^{\prime}})_{j}=(\mathbf{z}_{\alpha,\alpha^{\prime}})_{\beta}$, $(\mathbf{z}^{jj^{\prime}})_{i}=(\mathbf{z}^{\beta,\beta^{\prime}})_{\alpha}$. Then, it simple to check that $X_{ij}=X_{\alpha,\beta}$ satisfies the optimality conditions of equation A under conditions of the theorem and noticing that by the root-means-square and arithmetic mean (RMS-AM) inequality, we also have

[TABLE]

[TABLE]

∎

Lemma A.2**.**

Suppose that the ideal optimization in equation A has a solution where $X_{\alpha,\pi(\alpha)}>0$ holds for every $\alpha$. For every $\bm{\delta}=(\delta_{\alpha})_{\alpha}$ satisfying $\mathbf{1}^{T}\bm{\delta}=0$ and any choice of the optimal dual parameters $\{p_{\alpha},q_{\beta}\}$ we have that

[TABLE]

where $\delta^{\beta}=-\delta_{\pi^{-1}(\beta)}$. As a result in this case, equation A has optimal dual parameters $\{p_{\alpha},q_{\beta}\}$ satisfying

[TABLE]

Proof.

Denote the minimum value of $X_{\alpha,\pi(\alpha)}$ by $\epsilon$. Without loss of generality, we assume that $\|\bm{\delta}\|_{1}-\|\bm{\delta}\|_{\infty}\leq\epsilon$. Take $\alpha_{0}\in\arg\min\limits_{\alpha}|\delta_{\alpha}|$. Hence, $\|\bm{\delta}\|_{1}-\|\bm{\delta}\|_{\infty}=\sum\limits_{\alpha\neq\alpha_{0}}|\delta_{\alpha}|$.

Denote the optimal value of equation A by $C_{0}$. From the strong duality theorem we have that

[TABLE]

Take

[TABLE]

We notice that $\{p_{\alpha},q_{\beta}\}$ are feasible dual vectors for equation A.2. Hence, from the weak duality theorem we have

[TABLE]

[TABLE]

Now take the solution

[TABLE]

It is simple to check that $Y^{\prime}_{\alpha,\beta}$ is feasible in equation A.2. Moreover, we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We conclude that

[TABLE]

which proves the first part. Now, notice that for any pair $(\alpha_{1},\alpha_{2})$ of distinct indices, taking $\delta_{\alpha_{1}}=1$ and $\delta_{\alpha_{1}}=-1$ gives

[TABLE]

switching $\alpha_{1},\alpha_{2}$ yield

[TABLE]

Now, notice that from the optimality of equation A we have $p_{\alpha}+q_{\alpha}=D_{\alpha,\alpha}$, which leads to

[TABLE]

which yield

[TABLE]

The result is obtained by noticing that the set of optimal dual solutions is invariant under shift, i.e. $p_{i}+\lambda$ and $q_{i}-\lambda$ are also solutions for any $\lambda\in\mathbb{R}$. Hence, we may take $\lambda$ such that

[TABLE]

and hence

[TABLE]

Triangle inequality gives the result.

∎

A.1 Proof of Theorem 3.2

The first claim that $X_{ij}=X_{\alpha,\beta}$ follows by specializing part 2 of Theorem A.1 for the conditions of Theorem 3.2 with $a=1/2,b=c=0$: $\theta=\infty$ and inside clusters we have $R_{i,i^{\prime}}=S_{j,j^{\prime}=1}=1$, $R_{i,\alpha^{\prime}}=R_{i^{\prime},\alpha}$ and $S_{j,\beta^{\prime}}=R_{j^{\prime},\beta}$. Moreover $n_{\alpha}=m_{\beta}=m$.

Part 1 in the main text also is achieved by specializing part 1.a.: We will have $\sigma_{\alpha}=\omega_{\alpha}$, $n_{\alpha}=m_{\beta}=m$, $R_{\alpha,\alpha^{\prime}}=m^{2}R$ and $S_{\beta,\beta^{\prime}}=m^{2}$.

Finally, part 2 is a result of 1.b. with $\theta=\infty,\bm{\delta}=\mathbf{0}$.

A.2 Proof of Theorem 3.1

Based on theorem 3.2, we present a sketch of the proof for theorem 3.1. Under the assumptions of theorem 3.1, we directly verify that $\Lambda=\Lambda_{\alpha,\beta}=\nicefrac{{\sqrt{2}}}{{K(1+R)}}$. For part 1) $\delta=D-d-2\omega,\ \Delta\leq 2\omega$ and for part 2) $\delta=D^{2}-d^{2}$ are valid choices. Finally, for part 2) we may conclude by Chernoff bound that with a probability exceeding $1-\nicefrac{{1}}{{n^{10}}}$ (the power 10 is arbitrary) we have $\Delta=O(\omega\sqrt{(E+\omega)^{2}+1}\log(nK))$. Replacing these expression in the first part of theorem 3.2 gives us the result.

Appendix B Proof of Theorem 4.1

To simplify the notation, we introduce $\phi_{P+q}=I_{S_{q}}$ for $q=1,2\ldots,Q$, where $I_{S}$ denotes the indicator function of a convex set $S$. It is well-known that the proximal operator of $I_{S}$ coincides with the orthogonal projection operator onto $S$. Hence, we may simplify our algorithm to

[TABLE]

[TABLE]

where we introduce $\mathbf{g}_{P+q}=\mathbf{h}_{q}$ for $q=1,2\ldots,Q$ and denote $\bar{\mathbf{g}}_{t}=\sum\limits_{r=0}^{R}\mathbf{g}_{r}$ with $R=P+Q$. Moreover, $r_{t}$ is equal to either $p_{t}$ or $P+q_{t}$, depending on the random choice. We also define $\mathbf{x}^{\dagger}_{r,t}=\mathrm{prox}\left(\mathbf{x}_{t}+\mu\mathbf{g}_{r}\right)$ and hence $\mathbf{x}_{t+1}=\mathbf{x}^{\dagger}_{r,t}$.

To prove convergence, we adopt a so-called Lyaponov function approach. We introduce two non-negative functions $L,M$ of the state variables, $\mathbf{x}$ and $\{\mathbf{g}_{r}\}$ such that

[TABLE]

where $L_{t},M_{t}$ denote the values of $L,M$ at the variables of the $t^{\textrm{th}}$ iteration. Then, summing these inequalities up to an arbitrary time $t$ gives

[TABLE]

which by the non-negativity of $L$ implies

[TABLE]

In particular, we take

[TABLE]

where

[TABLE]

and $\mathbf{g}_{r}^{*}\in\partial\phi_{r}(\mathbf{x}^{*})$ for $r\in[R]$ satisfy the monotone inclusion problem in equation 4.6 at $\mathbf{x}^{*}$, i.e $\sum\limits_{r}\mathbf{g}_{r}^{*}=0$. Then, the non-negativity of each summand of $F_{t}$ follows from the convexity of $\phi_{r}$. The third term of $M_{t}$ is also positive for $\alpha<\frac{1+\sqrt{17}}{4}(1-\rho)$. This establishes the non-negativity of $M_{t}$. We further define

[TABLE]

and take

[TABLE]

B.1 Proof of Theorem Under equation B.1

Let us first prove the theorem assuming equation B.1 holds true for the given $L,M$. Then equation B.3 also holds true and we conclude from the definition of $M$ that

[TABLE]

Now from equation B.5 and the triangle inequality, we conclude that

[TABLE]

Further since the functions are $\beta-$Lipschitz, we observe that $\|\mathbf{g}_{r}^{*}\|\leq\beta$. Hence,

[TABLE]

Similarly,

[TABLE]

We conclude that

[TABLE]

and hence,

[TABLE]

Using Jensen’s inequality and equation LABEL:eq:middle_bounds, we obtain

[TABLE]

This proves the first bound in part 1 noting that for a suitable constant $c$ only depending on $\rho,\alpha$ and the constant in equation 4.7

[TABLE]

For the second bound in part 1 note that for $r=P+1,P+2,\ldots,P+Q$, we have $\mathrm{dist}(\mathbf{x}_{t},S_{r})\leq\|\mathbf{x}_{r,t}^{\dagger}-\mathbf{x}_{t}\|$, since by definition $\mathbf{x}_{r,t}^{\dagger}\in S_{r}$. We conclude that

[TABLE]

For part 2, note that

[TABLE]

We conclude from equation B.3 that

[TABLE]

which completes the proof.

B.2 Proof of equation B.1

It remains to prove equation B.1. We first state the following intermediate results that characterize the term $\mathbb{E}[L_{t}]-\mathbb{E}[L_{t-1}]$ in equation B.1. They can be proven by direct substitution and hence the proofs are neglected.

Lemma B.1**.**

The average dynamics of $G_{t}$ is given by:

[TABLE]

Lemma B.2**.**

The average dynamics of $\Gamma_{t}$ is given by:

[TABLE]

Lemma B.3**.**

The average dynamics of $D_{t}$ is given by:

[TABLE]

Now, we state a crucial inequality that connects the dynamics of $L$ to $M$:

Lemma B.4**.**

The following inequality holds at every time:

[TABLE]

Proof.

From the definition of a proximal operator for $r=1,2,\ldots,P$, we have $\frac{\mathbf{x}_{t}-\bar{\mathbf{x}}^{\dagger}_{r,t}}{\mu}+\mathbf{g}_{r}\in\partial\phi_{r}(\bar{\mathbf{x}}^{\dagger}_{r,t})$. hence

[TABLE]

Adding and subtracting the term $\langle\mathbf{g}^{*}_{r},\bar{\mathbf{x}}^{\dagger}_{r,t}-\mathbf{x}_{t}\rangle$ and summing over $r\in[P]$ gives

[TABLE]

Now, note that by the definition of a projection operator for $r=P+1,P+2,\ldots,R$ we have $\frac{\mathbf{x}_{t}-\bar{\mathbf{x}}^{\dagger}_{r,t}}{\mu}+\mathbf{g}_{r}$ is normal to $S_{r}$ at $\bar{\mathbf{x}}^{\dagger}_{r,t}$. Since $\mathbf{x}^{*}\in S_{r}$, we have

[TABLE]

Similarly, we obtain that

[TABLE]

Summing (44) and (45) over $r=P+1,P+2,\ldots,R$ and adding to (43) gives the desired result. ∎

B.2.1 Combining Bounds

To obtain equation B.1, we combine the inequalities in the above four lemmas in the following way. Respectively multiplying equation LABEL:eq:LD, equation LABEL:eq:LG and equation LABEL:eq:LGamma by $\frac{R}{2\mu}$, $\frac{1}{2\rho\mu\alpha}$ and $\frac{R\mu}{2\rho}$ and adding to equation B.21 and after straightforward calculations, we obtain:

[TABLE]

[TABLE]

[TABLE]

By invoking Jensen’s inequality, we have,

[TABLE]

which yields the desired result.

B.3 Proof of Theorem 4.2

The proximal operator of $\phi_{k,\bm{\zeta},\bm{\eta}}$ is defined as

[TABLE]

We introduce a change of variables by $\mathbf{u}=\nicefrac{{(\mathbf{x}+\mathbf{y})}}{{2}},\mathbf{v}=\nicefrac{{(\mathbf{x}-\mathbf{y})}}{{2}}$. First note that $\|a\|^{2}+\|b\|^{2}=(\|a+b\|^{2}+\|a-b\|^{2})/2$. Hence

[TABLE]

[TABLE]

Furthermore,

[TABLE]

Hence equation B.26 can be written as

[TABLE]

[TABLE]

This is separable over $\mathbf{u}$ and $\mathbf{v}$, and can be analytically solved. We get

[TABLE]

[TABLE]

The result is obtained by setting $\mathbf{x}=\mathbf{u}+\mathbf{v},\mathbf{y}=\mathbf{u}-\mathbf{v}$.

Appendix C Additional experiments

C.1 Impact of SON-Regularizer

We investigate the models on a simple dataset, shown in Fig. 2. We illustrate the behavior of each model with respect to two different values of its regularization parameter (low and high) respectively at the first and the second row (low: $\lambda_{1}=0.01,\lambda_{2}=0.0$, high: $\lambda_{1}=10,\lambda_{2}=5$). In the low setting, instead of $\lambda_{2}=0.0$ any other small value also yields consistent results. The source data, target data and transported source data are respectively shown by yellow, blue and red points. Each column of the sub-figures in Fig. 2 corresponds to a particular model performance respectively OT-l1l2, OT-lpl1, OT-Sinkhorn and OT-SON (our model). We observe that OT-SON yields stable and consistent results for different values of its parameters. Moreover, the data points transported by the proposed model are always informative providing a good representation of the underlying classes. Whereas, the other OT models are sensitive to the values of their regularization parameters and might thus transport the source data to somewhere in the middle of the actual target data, or away from the actual classes in the target domain.

We next study the interesting case where the source and target domains do not include the same number of classes, as shown in Fig. 3. In this experiment, we assume that the source data contains three classes, whereas the target domain has only two classes. Using the same color code as in Fig.2, we see in Fig. 3 the target classes and the transported source classes to the target domain are shown in yellow, blue and red respectively corresponding to OT-l1l2, OT-lpl1, OT-Sinkhorn and OT-SON. We again illustrate the behavior of each model w.r.t. two different values of its regularization parameter (low and high) respectively at the first and the second row. We observe that among all different models, only OT-SON with an appropriate parameter is able to identify that the source and the target domains have different number of classes, and subsequently, matches the corresponding classes correctly. It maps the superfluous class to a space between the two matched classes. However, the other models assign the superfluous class to the two other classes and do not distinguish the presence of such an extra class in the source domain. This observation is consistent with the assumptions made in (Courty et al.,, 2017). The unbalanced method in (Chizat et al.,, 2018) might be relevant but its use is unclear to us. In the last column of Fig. 3, the heat maps show the mapping cost among different source and target classes, and as well as the transport map obtained by our algorithm (OT-SON with a high regularization). We observe that the transport map respects the class structure.

C.2 Experiments on path-based data

In Fig. 4, we investigate the different OT-based domain adaptation models on a commonly-used synthetic dataset, wherein the three classes have diverse shapes and forms (Chang and Yeung,, 2008). In particular, we consider the case where one of the source classes is absent in the target domain. With the same number of classes in the source and target domains, the different models perform equally well. Fig. 4(a) shows a case where the source data (yellow points) and the target data (blue points), differ in the fact that the target data is missing the upper left Gaussian cloud of points appearing in the source data. Fig. 4(b) shows the two source and target datasets, as well as the transported data by our model (OT-SON). The transported data points are shown in red. We observe that our method avoids mapping the source data of the missing class to any of the present classes in the target domain. The points with white interior are those not assigned to any class in the target. This thus leads to a better prediction of the target data. In the table of Fig. 4(c), we compare the accuracy scores of different models on the target data, where our model yields the highest score.

C.3 Unsupervised domain adaptation

In all prior experiments, we have assumed that the class labels of the source data are available. This setup is consistent with the study in (Courty et al.,, 2017). We consequently evaluate in a side study the fully unsupervised setting, i.e., the case where no class label is available for the source or the target data. We consider the setting used in Fig. 3 with, this time, no given class labels. While the other methods fail for this task, the OT-SON with proper parameterization (i.e., the setting shown in the second row and the forth column) yields meaningful and consistent results. Fig. 5 shows the OT-SON results and the consistency of transport costs and transport maps computed by OT-SON.

C.4 Early stopping of the optimization

We study the early stopping of our optimization procedure. We use the data in Fig. 3 and investigate the results with different number of epochs. Here, we employ the OT-SON with proper parameterization, i.e., the results shown in the forth column and the second row for OT-SON in Fig. 3. In the experiments in Fig. 3 we performed the optimization with $20$ epochs. Here, we study early stopping, i.e., we study the quality of results if we stop after a smaller number of epochs. According to the results in Fig. 6, we observe that even after a small number of epochs, we obtain reliable and stable results that represent well the ultimate solution. Such a property is very important in practice, as it can significantly reduce the heavy computations. Fig. 7 illustrates the transport maps for different number of epochs. The different transport maps at different number of epochs are consistent with the transport cost shown in the last row of Fig. 7.

C.5 Diverse classes in the source

We next study the case where two of the three source classes have the same label, as shown in Fig. 8. In the source data (shown by yellow), the left and the middle data clouds have the same class labels. This example shows why the transport based on only the pairwise distances between the source and target data is insufficient. In Fig. 8, the left plot corresponds to $\lambda_{1}=\lambda_{2}=0$, the middle plot corresponds to $\lambda_{1}=10,\lambda_{2}=0.01$, and the right plot corresponds to $\lambda_{1}=100,\lambda_{2}=0.01$. We observe that the left plot (with $\lambda_{1}=\lambda_{2}=0$) fails to perform a proper transport of the source data. On the other hand, with incorporating our proposed regularization, the two different classes (even-though one of them is diverse) are properly transported to the target domain. We observe this kind of transfer in both of the middle ($\lambda_{1}=10,\lambda_{2}=0.01$) and right ($\lambda_{1}=100,\lambda_{2}=0.01$) plots.

C.6 Fewer classes in the source

In the experiments of Fig. 3, we studied the case where the number of source classes is larger the number of target classes. Here, we consider an opposite setting: we assume two classes in the source and three classes in the target, as illustrated in Fig. 9. The source, target and transported data points are respectively shown by yellow, blue, and red. We use the same setting and parameters as in Fig. 3, i.e., the first row corresponds to low regularization and the second row to high regularization (low regularization: $\lambda_{1}=0.01,\lambda_{2}=0.0$, high regularization: $\lambda_{1}=10,\lambda_{2}=5$). We observe that similar to the results in Fig. 3, only OT-SON with high regularization prevents splitting the source data among all the three target classes. The last row in Fig. 9 indicates the consistency between the mapping costs and transport map for this setting (for OT-SON with high regularization).

Appendix D Uniqueness

An elementary question concerning any optimization formulation, including the Kantorovich problem and its regularization in equation 2.2.1, is the uniqueness of their optimal solution, and a standard method for verifying uniqueness is to establish strong convexity of the objective function. Even though it is seen that the objective in (2.2.1) is not strongly convex, we are nevertheless able to identify conditions, under which the solution still remains unique. For this, we develop an alternative approach, which is not only useful in our framework, but may also be generically used in many similar problems including a wide range of linear programming (LP) relaxation problems, and for this reason it is first presented. Our approach is based on the following definition:

Definition D.1**.**

We call a (global) optimal solution $\mathbf{X}_{0}$ of a convex optimization problem

[TABLE]

where $\mathscr{F}(\ldotp)$ is a convex function and $\mathscr{S}$ is a convex set, a resistant optimal point if adding a linear perturbation term $\langle\tilde{\mathbf{D}},\mathbf{X}\rangle$ with sufficiently small coefficients in $\tilde{\mathbf{D}}$ to the objective leads to an arbitrarily small perturbation of the solution $\mathbf{X}_{0}$. In mathematical terms for any open neighborhood ${\mathcal{N}}$ of $\mathbf{X}_{0}$ there exists an open neighborhood ${\mathcal{M}}$ of $\tilde{\mathbf{D}}$ around $\tilde{\mathbf{D}}=\mathbf{0}$ such that

[TABLE]

Accordingly, we have the following result:

Theorem D.2**.**

A resistant optimal point of a convex optimization problem is its unique optimal point.

Proof.

Suppose that there exists a different optimal point $\mathbf{X}^{\prime}$. Take $\mathbf{D}_{0}=\frac{\mathbf{X}_{0}-\mathbf{X}^{\prime}}{\|\mathbf{X}_{0}-\mathbf{X}^{\prime}\|}$, $r=\|\mathbf{X}_{0}-\mathbf{X}^{\prime}\|$ and $\tilde{\mathbf{D}}=\epsilon\mathbf{D}_{0}$ for arbitrary $\epsilon>0$. Further, define ${\mathcal{N}}$ as the ball of radius $\delta=r/2$ centered at $\mathbf{X}_{0}$. Note that for each $\mathbf{Y}\in{\mathcal{N}}$ we have

[TABLE]

[TABLE]

Now, note that $\langle\tilde{\mathbf{D}},(\mathbf{Y}-\mathbf{X}_{0})+(\mathbf{X}_{0}-\mathbf{X}^{\prime})\rangle\geq-\delta\epsilon+r\epsilon>0,$ which establishes

[TABLE]

Hence, ${\mathcal{N}}\cap\arg\min\limits_{\mathbf{X}\in\mathscr{S}}\mathscr{F}(\mathbf{X})+\langle\tilde{\mathbf{D}},\mathbf{X}\rangle=\emptyset$ and since $\epsilon=\|\tilde{\mathbf{D}}\|$ is arbitrarily small, we conclude that $\mathbf{X}_{0}$ is not a resistant optimal point. This contradicts the assumption and shows that the solution is unique. ∎

Theorem D.2 is a general way to establish uniqueness. In fact, we can show that the strong convexity condition is a special case of this result:

Theorem D.3**.**

If $\mathscr{F}$ is continuous and strongly convex, then the global minimal point of $\mathscr{F}$ over a convex set $\mathscr{S}$ is resistant.

Proof.

Denote the optimal point by $\mathbf{X}^{*}$. By strong convexity, there exists a $\gamma>0$ such that for any feasible point $\mathbf{X}\in\mathscr{S}$, we have $\mathscr{F}(\mathbf{X})-\mathscr{F}(\mathbf{X}^{*})\geq\frac{\gamma}{2}\|\mathbf{X}-\mathbf{X}^{*}\|_{\mathrm{F}}^{2}$. Take $\mathscr{G}=\mathscr{F}+\langle\tilde{\mathbf{D}},\mathbf{X}\rangle$ and note that $\mathscr{G}(\mathbf{X})-\mathscr{G}(\mathbf{X}^{*})\geq\frac{\gamma}{2}\|\mathbf{X}-\mathbf{X}^{*}\|_{\mathrm{F}}^{2}+\langle\tilde{\mathbf{D}},\mathbf{X}-\mathbf{X}^{*}\rangle\geq\frac{\gamma}{4}\|\mathbf{X}-\mathbf{X}^{*}\|_{\mathrm{F}}^{2}-\frac{2}{\gamma}\|\tilde{\mathbf{D}}\|_{\mathrm{F}}^{2}$. This shows that $\mathscr{G}>\mathscr{G}(\mathbf{X}^{*})$ and hence does not have any global optimal point outside the closed sphere $\{\mathbf{X}\mid\|\mathbf{X}-\mathbf{X}^{*}\|_{\mathrm{F}}\leq\frac{\sqrt{8}}{\gamma}\|\tilde{\mathbf{D}}\|_{\mathrm{F}}\}$. Since $\mathscr{G}$ is continuous, it also attains a minimum inside the sphere, which then becomes the global optimal point. We conclude that for any $\epsilon>0$, taking $\|\tilde{\mathbf{D}}\|<\frac{\gamma\epsilon}{\sqrt{8}}$ leads to an optimal solution inside a ball of radius $\epsilon$ centered at $\mathbf{X}^{*}$. This shows that the solution is resistant. ∎

Uniqueness for equation 2.2.1: One special case of resistant optimal points, that will be useful in our analysis, is when there exists a neighborhood ${\mathcal{M}}$ of $\mathbf{0}$ such that

[TABLE]

We call such a resistant optimal point an extremal optimal point. Later, we consider an analysis where we give conditions on $\mathbf{D}$ to ensure that a desired solution $\mathbf{X}^{*}$ is achieved. Our strategy for uniqueness in this analysis is to show that under the same conditions, the desired optimal point is also extremal and hence unique, according to Theorem 1. In the case of the problem in equation 2.2.1, adding the term $\langle\tilde{\mathbf{D}},\mathbf{X}\rangle$ modifies the cost matrix $\mathbf{D}$ to $\mathbf{D}+\tilde{\mathbf{D}}$. Hence, being an extremal optimal point is in this case equivalent to the solution $\mathbf{X}^{*}$ being maintained following a perturbation of the matrix $\mathbf{D}$ in a sufficiently small open neighborhood. This is easy to achieve in our planted model analysis, because the optimality of $\mathbf{X}^{*}$ is guaranteed by a set of inequalities on $\mathbf{D}$, which remain valid under small perturbations, simply by requiring the inequalities to be strict. As seen, Theorem D.2 and extremal optimality, in particular, can be powerful tools for establishing uniqueness beyond strong convexity.

Appendix E Remarks on the Optimization Algorithm

Efficient Computation: While the objective in (LABEL:eq:fullobj) may appear complex as it involves $n^{2}$ terms, the associated algorithm is stochastic and incremental, thus only involving one term in (LABEL:eq:fullobj) for each iteration, thus greatly reducing the complexity as a result. The simplification of the algorithm is also due to the proximal update detailed in Theorem 4.2 (and subsequent projection) used in each iteration update of a pair of rows or columns. We further note that an early stopping typical of stochastic schemes is likely, making a full-run to convergence unnecessary (see Section 4.4 in (Bottou et al.,, 2018)), and in practice avoiding the impact of the $n^{2}$ terms on the performance. When the underlying data satisfies the structure of the stochastic block model, the problem size is essentially $B^{2}\ll n^{2}$, as the number of required iterations is determined by an adequate sampling of all blocks.

Just-in-Time Update: In our problem of interest in equation 2.2.1, the number of variables quadratically grows with the problem size. For such problems, incremental algorithms may become infeasible in large-scale. Note that each iteration of our algorithm includes proximal and projection operators, that update only a small group of variables. This allows us to apply the Just-in-Time approach in (Schmidt et al.,, 2017) to resolve the problem with the number of variables. In our problem, each term $\phi_{n}(x)$ and constraint $S_{m}$ only involves a small subset $x_{I_{n}}:=(x_{i},i\in I_{n})$ of the variables, where $I_{n}\subseteq[D]$. Hence, the projection and proximal operators alter only a small subset of variables, dramatically reducing the amount of computation. We exploit this to give an algorithm that has much cheaper per-iteration cost. Note that the vanilla algorithm explained in equation 4.4 and equation 4.5 still operates on the full set of variables as the memory vectors become non-sparse by the updating rule in equation 4.5. We resolve this issue by following the Just-in-Time approach in (Schmidt et al.,, 2017) and modifying equation 4.5 to

[TABLE]

where $I_{t}$ denotes the set of variables involved in the $t^{\textrm{th}}$ iteration and we define $(\mathbf{y})_{I}$ for a vector $\mathbf{y}=(y_{1},y_{2},\ldots,y_{d})$ as a vector $\mathbf{y}^{\prime}=(y_{1}^{\prime},y_{2}^{\prime},\ldots,y_{d}^{\prime})$ such that

[TABLE]

where $K=M+N$ and $K_{i}$ is the number of objective terms $\phi_{n}$ and constraint sets $S_{m}$ including the $i^{\textrm{th}}$ variable $x_{i}$.