\subsubsection*{\bibname}}

\usepackage{hyperref}
\usepackage{url}

\usepackage{xcolor} \usepackage{graphicx} \usepackage{booktabs} \usepackage{array} \usepackage{mathtools} \usepackage{amsthm, amssymb, amscd, amsfonts} \usepackage[ampersand]{easylist} \usepackage{multirow} \usepackage{booktabs} \usepackage{bm} \usepackage{bbm} \usepackage{dcolumn} \usepackage{caption}
\usepackage{subcaption} \usepackage{epstopdf} \usepackage{comment} \usepackage[export]{adjustbox}
\usepackage{hyperref}
\usepackage{wrapfig}
\usepackage{thm-restate} \usepackage{xspace} \usepackage[normalem]{ulem}

\usepackage{enumitem}
\setitemize{noitemsep,topsep=0pt,parsep=0pt,partopsep=0pt}
\setenumerate{noitemsep,topsep=0pt,parsep=0pt,partopsep=0pt}
\usepackage{hyperref}
\usepackage{algorithm} \usepackage{algorithmic} %
\newtheorem{theorem}{Theorem}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{definition}{Definition}
\newtheorem{intuition}{Intuition}
\theoremstyle{definition}
\newtheorem{observation}{Observation}
\theoremstyle{definition}
\newtheorem{hypothesis}{Hypothesis}
\theoremstyle{definition}
\newtheorem{remark}{Remark}
\theoremstyle{definition}
\newtheorem{conjecture}{Conjecture}
\theoremstyle{definition}
\newtheorem{question}{Question}

\renewcommand{\definitionautorefname}{Def.}
\renewcommand{\propositionautorefname}{Proposition}
\renewcommand{\remarkautorefname}{Remark}
\renewcommand{\lemmaautorefname}{Lemma}
\renewcommand{\corollaryautorefname}{Cor.}
\renewcommand{\equationautorefname}{Eqn.}
\renewcommand{\figureautorefname}{Fig.}
\newcommand{Alg.}{Alg.}

\renewcommand{\algorithmicrequire}{\textbf{Input:}}
\renewcommand{\algorithmicensure}{\textbf{Output:}}

\DeclareMathOperator*{\argmin}{arg\,min\,}
\DeclareMathOperator*{\argmax}{arg\,max\,}
\DeclareMathOperator*{\Id}{id}
\DeclareMathOperator*{\proj}{proj}
\DeclareMathOperator*{\linspan}{span}

\newcommand{p}{p}
\newcommand{\mathbb{E}}{\mathbb{E}}
\newcommand{\,|\,}{\,|\,}
\newcommand{, }{, }

\newcommand{\ensuremath{\mathbb{Z}}}{\ensuremath{\mathbb{Z}}}
\newcommand{\ensuremath{\Z_{+}}}{\ensuremath{\ensuremath{\mathbb{Z}}_{+}}}
\newcommand{\ensuremath{\Z_{-}}}{\ensuremath{\ensuremath{\mathbb{Z}}_{-}}}
\newcommand{\ensuremath{\mathbb{R}}}{\ensuremath{\mathbb{R}}}
\newcommand{\ensuremath{\R_{+}}}{\ensuremath{\ensuremath{\mathbb{R}}_{+}}}
\newcommand{\ensuremath{\R_{-}}}{\ensuremath{\ensuremath{\mathbb{R}}_{-}}}

\newcommand{\mathcal{L}}{\mathcal{L}}

\newcommand{n}{n}
\newcommand{\nobs}{n}
\newcommand{d}{d}
\newcommand{\nvar}{d}
\newcommand{m}{m}
\newcommand{\nlayer}{m}
\newcommand{\nlayers}{\nlayer}
\newcommand{J}{J}
\newcommand{k}{k}

\newcommand{\ij}{\ij}
\newcommand{\it}{\it}
\newcommand{i}{i}
\renewcommand{\ij}{j}
\newcommand{s}{s}
\renewcommand{\it}{t}
\newcommand{\ell}{\ell}
\newcommand{k}{k}

\newcommand{\u}{\u}

\newcommand{x}{x}
\newcommand{\bm{\x}}{\bm{x}}
\newcommand{X}{X}
\newcommand{\mathcal{\xmat}}{\mathcal{X}}
\newcommand{y}{y}
\newcommand{\bm{\y}}{\bm{y}}
\newcommand{Y}{Y}
\newcommand{\mathcal{\ymat}}{\mathcal{Y}}
\newcommand{z}{z}
\newcommand{\bm{\z}}{\bm{z}}
\newcommand{Z}{Z}
\newcommand{\mathcal{\zmat}}{\mathcal{Z}}

\newcommand{\bm{a}}{\bm{a}}
\newcommand{A}{A}
\newcommand{\mathcal{\amat}}{\mathcal{A}}
\newcommand{\mathcal{\amat}^{\complement}}{\mathcal{A}^{\complement}}

\newcommand{\bm{b}}{\bm{b}}
\newcommand{B}{B}
\newcommand{\mathcal{\bmat}}{\mathcal{B}}
\newcommand{\mathcal{\bmat}^{\complement}}{\mathcal{B}^{\complement}}

\newcommand{\widehat{x}}{\widehat{x}}
\newcommand{\widehat{\bm{x}}}{\widehat{\bm{x}}}
\newcommand{\widehat{X}}{\widehat{X}}
\newcommand{\widehat{y}}{\widehat{y}}
\newcommand{\widehat{\bm{y}}}{\widehat{\bm{y}}}
\newcommand{\widehat{Y}}{\widehat{Y}}
\newcommand{\widehat{z}}{\widehat{z}}
\newcommand{\widehat{\bm{z}}}{\widehat{\bm{z}}}
\newcommand{\widehat{Z}}{\widehat{Z}}
\newcommand{\widehat{u}}{\widehat{u}}
\newcommand{\widehat{\bm{u}}}{\widehat{\bm{u}}}
\newcommand{\widehat{U}}{\widehat{U}}
\newcommand{\widehat{v}}{\widehat{v}}
\newcommand{\widehat{\bm{v}}}{\widehat{\bm{v}}}
\newcommand{\widehat{V}}{\widehat{V}}

\renewcommand{\u}{u}
\newcommand{\bm{\u}}{\bm{\u}}
\newcommand{U}{U}
\newcommand{\mathcal{\umat}}{\mathcal{U}}
\renewcommand{\v}{v}
\newcommand{\bm{\v}}{\bm{\v}}
\newcommand{V}{V}
\newcommand{\mathcal{\vmat}}{\mathcal{V}}
\renewcommand{p}{p}
\newcommand{\bm{\prob}}{\bm{p}}
\newcommand{P}{P}
\newcommand{\mathcal{\probmat}}{\mathcal{P}}

\newcommand{w}{w}
\newcommand{\bm{\w}}{\bm{w}}
\newcommand{W}{W}
\newcommand{\mathcal{\wmat}}{\mathcal{W}}

\newcommand{\mathcal{C}}{\mathcal{C}}
\newcommand{\pi}{\pi}
\newcommand{\bm{\mu}}{\bm{\mu}}
\newcommand{\mathcal{N}}{\mathcal{N}}
\newcommand{\mathbf{I}}{\mathbf{I}}

\newcommand{\mathbb{P}}{\mathbb{P}}
\newcommand{{\hat{p}}}{{\hat{p}}}
\newcommand{{\mathcal{M}}}{{\mathcal{M}}}

\newcommand{\pi}{\pi}
\newcommand{\bm{\tq}}{\bm{\pi}}
\newcommand{\Pi}{\Pi}
\newcommand{\mathcal{\tqmat}}{\mathcal{\Pi}}

\newcommand{\tilde{q}}{\tilde{q}}
\newcommand{\tilde{\bm{q}}}{\tilde{\bm{q}}}
\newcommand{\hat{Q}}{\hat{Q}}
\newcommand{\hat{\mathcal{Q}}}{\hat{\mathcal{Q}}}

\newcommand{\textnormal{KL}}{\textnormal{KL}}
\newcommand{\textnormal{diag}}{\textnormal{diag}}
\newcommand{\text{T}} \newcommand{\supp}{\text{supp}}{\text{T}} \newcommand{\supp}{\text{supp}}
\newcommand{\textrm{d}} \newcommand{\E}{\mathbb{E}}{\textrm{d}} \newcommand{\E}{\mathbb{E}}
\newcommand{\mathcal{F}} \newcommand{\wass}{\mathcal{W}}{\mathcal{F}} \newcommand{\wass}{\mathcal{W}}

\newcommand{\Psi}{\Psi}
\newcommand{\widetilde{\deepfunc}}{\widetilde{\Psi}}
\newcommand{\mathcal{I}}{\mathcal{I}}
\newcommand{\R}{\ensuremath{\mathbb{R}}}
\newcommand{\omega}{\omega}
\newcommand{\Omega}{\Omega}
\newcommand{Q}{Q}
\newcommand{q}{q}
\newcommand{F}{F}
\newcommand{\bm{F}}{\bm{F}}
\newcommand{k}{k}
\newcommand{\mathcal{I}}{\mathcal{I}}
\newcommand{\mathcal{H}}{\mathcal{H}}
\newcommand{D_{\text{KL}}}{D_{\text{KL}}}
\newcommand{m}{m}
\newcommand{b}{b}
\newcommand{k}{k}

\newcommand\independent{\protect\mathpalette{\protect\independenT}{\perp}}
\def\independenT#1#2{\mathrel{\rlap{$#1#2$}\mkern2mu{#1#2}}}

\newcommand{R}{R}
\newcommand{f}{f}
\newcommand{\bm{f}}{\bm{f}}
\newcommand{\mathcal{F}}{\mathcal{F}}
\newcommand{g}{g}
\newcommand{\bm{g}}{\bm{g}}
\newcommand{\mathcal{G}}{\mathcal{G}}

\newcommand{\textnormal{JSD}}{\textnormal{JSD}}
\newcommand{G}{G}
\newcommand{D}{D}
\newcommand{\D}{D}
\newcommand{\mathcal{L}}{\mathcal{L}}
\newcommand{\mathcal{I}}{\mathcal{I}}

\newcommand{\hspace{-0.5em}}{\hspace{-0.5em}}
\newcommand{\tabeqn}[1]{\( \hspace{-0.5em} \begin{array}[t]{l} \displaystyle #1 \end{array} \)}

\newcommand{d}{d}
\newcommand{\mathcal{S}}{\mathcal{S}}

\newcommand{\mathbbm{1}}{\mathbbm{1}}

\renewcommand{d}{d}
\renewcommand{\mathcal{L}}{\mathcal{L}}

\newcommand{H_{\textnormal{C}}}{H_{\textnormal{C}}}
\newcommand{\hl}[1]{{\color{red} #1}}

\newcommand{\mathcal{P}}{\mathcal{P}}
\newcommand{\mathcal{D}}{\mathcal{D}}
\newcommand{\tagcomment}[1]{\\ & \quad\quad \text{(#1) \notag}}
\newcommand{R}{R}
\newcommand{\tilde{R}}{\tilde{R}}

\newcommand{\mathcal{N}(0,1)}{\mathcal{N}(0,1)}
\newcommand{\textnormal{bary}}{\textnormal{bary}}
\newcommand{{\textstyle \frac{1}{2}}}{{\textstyle \frac{1}{2}}}
\newcommand{1/2}{1/2}

\renewcommand{\text{T}} \newcommand{\supp}{\text{supp}}{T}
\newcommand{t}{t}
\renewcommand{D}{D}
\newcommand{d}{d}
\newcommand{\textnormal{pre}}{\textnormal{pre}}
\newcommand{\textnormal{post}}{\textnormal{post}}
\newcommand{\textnormal{mid}}{\textnormal{mid}}
\newcommand{\W}{W}

\newcommand{\textnormal{MaxEnt}}{\textnormal{MaxEnt}}
\newcommand{\textnormal{Ind}}{\textnormal{Ind}}
\newcommand{\textnormal{DensDestructor}}{\textnormal{DensDestructor}}

\newcommand{\textnormal{SIE}\xspace}{\textnormal{SIE}\xspace}
\newcommand{k}{k}
\newcommand{\textnormal{na\"ive}\xspace}{\textnormal{na\"ive}\xspace}
\newcommand{\textnormal{Na\"ive}\xspace}{\textnormal{Na\"ive}\xspace}

\newcommand{\wrapproof}[1]{
	\iftrue #1
	\fi }

\newcommand{w}{w}
\newcommand{\bm{\weight}}{\bm{w}}
\renewcommand{d}{d}
\newcommand{X}{X}
\newcommand{\textnormal{tree}}{\textnormal{tree}}
\newcommand{\textnormal{Pa}}{\textnormal{Pa}}
\newcommand{\mathcal{L}}{\mathcal{L}}
\newcommand{\mathcal{R}}{\mathcal{R}}
\newcommand{s}{s}
\newcommand{t}{t}
\renewcommand{p}{p}
\newcommand{m}{m}
\newcommand{\textnormal{hist}}{\textnormal{hist}}

\newcommand{\textbf{return}}{\textbf{return}}

\newcommand{\todo}[1]{{\color{red}TODO: #1}}
\newcommand{\david}[1]{{\color{blue}DI: #1}}
\newcommand{\add}[1]{{\color{magenta}Add: #1}}
\newcommand{\edit}[1]{{\color{red}EDIT: #1}}
\newcommand{\doubts}[1]{{\color{red}DOUBTS: #1}}
\newcommand{\zeyu}[1]{{\color{green}#1}}
\newcommand{\mysubsubsection}[1]{\paragraph{#1.}}
\renewcommand{\todo}[1]{}
\renewcommand{\david}[1]{}

\title{Iterative Barycenter Flows}

\author{
	David I. Inouye\\
	Purdue University\\
	\texttt{dinouye@purdue.edu} \\
	\And
	Zeyu Zhou\\
	Purdue University\\ \texttt{zhou1059@purdue.edu} \\
	\And
	Ziyu Gong\\
	Purdue University\\
	\texttt{gong123@purdue.edu}\\
	\And
	Pradeep Ravikumar\\
	Carnegie Mellon University\\
	pradeepr@cs.cmu.edu
}
\begin{document}
	\maketitle

	\begin{abstract}
		The task of mapping two or more distributions to a shared representation has many applications including fair representations, batch effect mitigation, and unsupervised domain adaptation.
		However, most existing formulations only consider the setting of two distributions, and moreover, do not have an identifiable, unique shared latent representation.
		We use optimal transport theory to consider a natural multiple distribution extension of the Monge assignment problem we call the symmetric Monge map problem and show that it is equivalent to the Wasserstein barycenter problem.
		Yet, the maps to the barycenter are challenging to estimate.
		Prior methods often ignore transportation cost, rely on adversarial methods, or only work for discrete distributions.
		Therefore, our goal is to estimate invertible maps between two or more distributions and their corresponding barycenter via a simple iterative flow method.
		Our method decouples each iteration into two subproblems: 1) estimate simple distributions and 2)
		estimate the invertible maps to the barycenter via known closed-form OT results.
		Our empirical results give evidence that this iterative algorithm approximates the maps to the barycenter.
	\end{abstract}

	\section{INTRODUCTION}
	The task of mapping two or more distributions to a shared representation has attracted increasing interest due to its varied applications. These include fair representations \citep{zemel2013learningfair}, batch effect mitigation \citep{haghverdi2018batch}, unsupervised domain adaptation \citep{hu2018duplex}, and generative models \citep{grover2020alignflow}.
	For example, \citet{zemel2013learningfair} estimate a shared latent representation of the class-conditional distributions that simultaneously obfuscates any information about protected attributes (e.g., race) while preserving all other information useful for classification.
	For genetic data, \citet{haghverdi2018batch} attempt to mitigate batch effects (i.e., irrelevant shifts in the data between batches caused by non-biological factors) by estimating a shared representation among batches;
	this enables the integration and analysis of multiple datasets collected at different laboratories.
	\citet{hu2018duplex} perform unsupervised domain adaptation by mapping the source and target domains to a shared latent representation.
	In most existing formulations for this task (of finding mappings from multiple distributions to a shared representation) however, the mapping to the shared representation and sometimes even the shared representation itself is under constrained and non-unique (e.g., \citep{grover2020alignflow}). This lack of identifiability of course has algorithmic implications (such as possible oscillatory non-convergent behavior), but moreover also poses theoretical concerns (since the target of the estimation procedures is unclear). While it might be possible to impose constraints on the mappings to ensure identifiability, it is not clear what such natural constraints might be.

	In this paper, we propose to use the theory of optimal transport (OT) towards the task of \emph{identifiable mappings} from multiple distributions to a shared representation: by adding the additional restriction of requiring the mappings to have the \emph{least transportation cost}. This can be shown to provide a natural and importantly --- unique --- shared latent representation.
	Moreover, as we show in the experiments, the resulting shared latent representations are much more meaningful: because they aim to minimize transportation cost, the maps are parsimonious in transforming each distribution, and in particular, preserve shared structure among the different distributions.
	For instance, the shared representation among digits should be digit-like. Whereas, methods ignore transportation cost, typically transform the class distributions to some baseline distribution such as uniform or Gaussian.
	We believe that such un-parsimonious transformations would accordingly incur greater sample complexity;
	this is particularly so if the distributions are approximately singular and lie on a lower dimensional manifold (e.g., images of digits).

	Recall the classical Monge assignment problem which aims to find the mapping from one distribution to another with the least cost. One (minor) caveat with this formulation is that this is not a symmetric formulation, but more importantly, this is focused on the case of two distributions. We thus consider a natural multiple distribution symmetric extension of the Monge assignment problem, which we call the symmetric Monge map problem: where we aim to find mappings from multiple distributions to a shared latent distribution, but with the least transport costs.  Interestingly, under some regularity conditions, this turns out to be equivalent to the \emph{barycenter} of the given distributions, which can be seen as an extension of the Euclidean mean to distributions \citep{peyre2019computational}.
	From OT theory, we have that the optimal map from each of the distributions to this barycenter for continuous distributions is unique and invertible. Thus, OT theory identifies a principled and theoretically-grounded shared representation with unique corresponding maps.

	There is however a critical computational caveat with the formulation above: barycenter distributions and OT maps are challenging to estimate, especially given only samples.
	Most prior practical barycenter work has focused on estimating barycenters in discrete spaces (see e.g., \citep{peyre2019computational})
	While the continuous problem can be approximated by discretizing the space via histograms or clustering, this approach will fail in higher dimensions because of the curse of dimensionality.
	Another approach is to solve the barycenter problem directly on samples and then use a neural network to interpolate for new data points \citep{seguy2018large}, but this process only learns the map indirectly and does not guarantee invertibility.
	Thus, in spite of the theoretical advantage of identifiability and the parsimony of shared representations, simply using existing Barycenter algorithms might be less scalable than the non-OT deep generative approaches that are popularly used for this task.
	In particular, flow-based deep generative methods learn invertible mappings between distributions by construction (e.g., \citep{dinh2017density,kingma2016improved,kingma2018glow}), but these methods do not learn an identifiable mapping, ignore transportation costs (e.g., \citep{dinh2017density}), and moreover, might require adversarial methods (e.g., \citep{grover2020alignflow}), which can be difficult to optimize in practice \citep{lucic2018gans,kurach2019the}.

	We thus aim to develop a simple class of scalable iterative algorithms for approximating the maps to the shared barycenter distribution.
	Instead of trying to solve the global problem directly, we propose to iteratively solve simpler subproblems by first approximating the class distributions and then estimating the barycenter maps by leveraging known closed-form results (particularly for Gaussians and 1D problems).
	This can be seen as optimizing approximate subproblems rather than approximately solving a large joint problem---similar in spirit to Newton's method or coordinate descent.
	We summarize our contributions as follows:
	%
	\begin{itemize}
		\item We formulate the problem of estimating the maps from multiple distributions to a shared latent representation as a symmetric Monge map problem whose solution is unique; we show that it is equivalent to the barycenter problem and uncover its relationship to prior adversarial and flow-based methods.
		\item We propose a general iterative meta-algorithm to approximately optimize this symmetric Monge map problem via iteratively solving two steps:
		given respective samples, estimate distributions from simple parametric families, estimate the symmetric Monge map problem for these estimated distributions via known closed-form results, and then using the Monge maps to transform the respective samples for the next iteration.
		\item We provide three concrete instantiations of the above general meta-algorithm and demonstrate that our approach approximates the mappings to the barycenter on several simulated and real-world datasets.
	\end{itemize}

	\section{SYMMETRIC MONGE MAPS}

	\subsection{Background: Classical Monge Maps}
	We will first review some standard OT definitions. First, we note that $T_{\sharp} \alpha$ will denote the push-forward measure of $\alpha$ via the map $T$ (see \citep[Def.~2.1]{peyre2019computational}).
	\begin{definition}[Monge problem {\citep[Remark~2.7]{peyre2019computational}} ]
		\label{def:monge-problem}
		Given two measures $(\alpha, \beta)$ supported on two spaces $(\mathcal{\xmat}, \mathcal{\ymat})$ and a cost function $c(\cdot, \cdot)$, the Monge problem is defined as finding the map $T \colon \mathcal{\xmat} \to \mathcal{\ymat}$ that solves:
		\begin{align*}
		\begin{aligned}
		\arg\min_{T} &\int_{\mathcal{\xmat}} c(x, T(x)) d\alpha(x) \quad \textnormal{s.t.} \quad T_{\sharp} \alpha = \beta \, .
		\end{aligned}
		\end{align*}
	\end{definition}
	We will denote the optimal Monge map between two distributions as $T^*_{\alpha \to \beta}$ where $\alpha$ and $\beta$ represent two distributions.

	\subsection{Symmetric Monge Map Formulation}
	We note that the Monge formulation is asymmetric because the two distributions have distinct roles.
	While in theory the role of the distributions does not matter because $T^*_{2\to1} \equiv (T^{*}_{1\to2})^{-1}$, in practice the estimated map $\hat{T}$ will vary depending on which distribution is the source distribution.
	Finally and more importantly, the Monge problem in its original formulation only considers two distributions but we want to consider more than two distributions.
	\begin{definition}[Symmetric Monge Map]
		\label{def:symmetric-monge}
		Given a set of continuous input measures $(\mu_1, \dots, \mu_k)$ defined on some continuous space $\mathcal{\xmat}$, a non-negative weight vector $\bm{\weight} \geq 0$ such that $\sum_j w_j = 1$, and cost function $c(\cdot, \cdot)$, the \emph{symmetric Monge map problem} is defined as:
		\begin{align}
		\argmin_{T_1, T_2, \cdots, T_k} & \quad \sum_{j=1}^k w_j \int_{\mathcal{\xmat}_j} c(x, T_j(x)) d\mu_j(x) \label{eqn:transportation-cost} \\
		\textnormal{s.t.} &\quad T_{j\sharp} \mu_j = T_{j'\sharp} \mu_{j'} \,\,\, \forall j \neq j'  \, ,
		\label{eqn:pushforward-condition}
		\end{align}
		where \autoref{eqn:transportation-cost} is the transportation cost and \autoref{eqn:pushforward-condition} is known as the \emph{pushforward condition}.
	\end{definition}
	The original Monge problem can be recovered if $T_2 = \Id$ and $w_2=0$.  Thus, this problem can be seen as a symmetric relaxation of the Monge problem for two or more distributions.
	The pushforward condition in \autoref{eqn:pushforward-condition} is hard to satisfy for continuous distributions known only through samples.

	\subsection{Equivalence to Barycenter}
	We prove that our symmetric Monge map problem is equivalent to the barycenter problem (proof in appendix).
	We first review the definitions of the Kantorovich relaxation and the barycenter distribution, and then present our theorem.
	For this paper, we will assume $c(x,y) = \|x-y\|_2^2$ and that one of the measures has a density so that the barycenter is unique \citep{agueh2011barycenters}.
	\begin{definition}[Kantorovich Relaxation {\citep[Remark~2.13]{peyre2019computational}}]
		\label{def:kantorovich-relaxation}
		Given the same variables as \autoref{def:monge-problem}, the Kantorovich problem is defined as:
		$    \mathcal{L}_c(\alpha, \beta) \triangleq \min_{\pi \in \mathcal{U}(\alpha, \beta)} \int_{\mathcal{\xmat} \times \mathcal{\ymat}} c(x, y) d\pi(x,y) \, ,$
		where $\pi$ is a joint distribution over $\mathcal{\xmat}$ and $\mathcal{\ymat}$ such that the marginals are equal to $\alpha$ and $\beta$ respectively (denoted by $\mathcal{U}(\alpha, \beta)$).
	\end{definition}

	\begin{definition}[Wasserstein Barycenter {\citep[Remark 9.1]{peyre2019computational}}]
		\label{def:wasserstein-barycenter}
		Given a set of input measures $(\mu_1, \cdots, \mu_k)$ defined on some space $\mathcal{\xmat}$, weights $\bm{\weight}$ such that $\sum_j w_j = 1$, the barycenter is defined as:
		$    \textnormal{bary}(\mu_1, \mu_2, \cdots, \mu_k; \bm{\weight}) \triangleq \argmin_{\nu} \sum_{j=1}^k w_j \mathcal{L}_c(\nu, \mu_j) \, ,$
		where $\mathcal{L}_c$ is defined in \autoref{def:kantorovich-relaxation}.
	\end{definition}
	\begin{theorem}[Barycenter Equivalence for Symmetric Monge Map]
		\label{thm:barycenter-equivalence}
		For $c(x,y) = \|x-y\|_2^2$ where the measures have densities, the symmetric Monge map solution (\autoref{def:symmetric-monge}) gives the optimal maps between the class distributions and the barycenter distribution (\autoref{def:wasserstein-barycenter}), i.e., $T^*_j = T^*_{j \to \nu}$ where $\nu = \textnormal{bary}(\mu_1, \mu_2, \dots, \mu_k; \bm{\weight})$.
	\end{theorem}

	\subsection{Relation to Prior Methods}

	\paragraph{Adversarial Methods}
	CycleGAN \citep{CycleGAN2017} minimizes the fairly complex objective function:
	\begin{align}
	&\argmin_{G,F} d_{\textnormal{adv}}(G_\sharp \mu_1, \mu_2) + d_{\textnormal{adv}}(F_\sharp \mu_2, \mu_1) \label{eqn:cycle-adv-terms}\\
	&\!+ \!\lambda_1 \!\Big(\E_{\mu_1}[\|F(G(x))\! -\! x\|_1] \!+\! \E_{\mu_2}[\|G(F(x)) \!-\! x\|_1]\!\Big) \label{eqn:cycle-consistency-terms}\\
	&\!+\! \lambda_2 \Big(\E_{\mu_1}[\|F(x) - x\|_1] + \E_{\mu_2}[\|G(x) - x\|_1]\Big) \,, \label{eqn:cycle-identity-terms}
	\end{align}
	where the distance $d_{\textnormal{adv}}$ \eqref{eqn:cycle-adv-terms} approximates the distance between distributions via adversarial learning (i.e., minimax learning) and approximates the pushforward constraint;
	the cycle consistency term \eqref{eqn:cycle-consistency-terms} can be seen as a relaxation of an invertibility constraint; and the ``identity'' regularization \eqref{eqn:cycle-identity-terms} can be seen as a transportation cost term with $c(x,y) = \|x-y\|_1$.
	Thus, CycleGAN is approximately solving the Monge problem.
	StarGAN \citep{choi2018stargan} generalizes CycleGAN to more than two domains.
	However, these approaches cannot guarantee invertibility and require expensive and challenging adversarial learning \citep{lucic2018gans,kurach2019the}.

	\paragraph{Flow Methods}
	AlignFlow \citep{grover2020alignflow} extends CycleGAN by using invertible models so that the cycle consistency constraint is satisfied by construction:
	\begin{align}
	\argmin_{T_1, T_2} &d_{\textnormal{adv}}((T^{-1}_{2} \!\circ\! T_{1})_\sharp \mu_1, \mu_2\!) \!+\! d_{\textnormal{adv}}((T^{-1}_{1} \circ T_{2})_\sharp \mu_2, \mu_1\!)
	+ \lambda \Big(\textnormal{KL}(\mu_1, T_{1\sharp}^{-1} \alpha) + \textnormal{KL}(\mu_2, T_{2\sharp}^{-1} \alpha)\Big) \,,
	\end{align}
	where the first two distance terms (equivalent to CycleGAN) are implemented using adversarial learning, $\alpha$ is a Gaussian prior distribution, and the KL terms are implemented via maximum likelihood.
	The first two terms are equivalent to $d(T_{j\sharp} \mu_j, T_{j'\sharp} \mu_{j'})$---which aligns with the symmetric Monge pushforward constraint \autoref{eqn:pushforward-condition}---if the distance is invariant under invertible transformations (e.g., JSD as approximated by adversarial learning $d_{\textnormal{adv}}$).
	However, unlike the symmetric Monge formulation, AlignFlow ignores transportation costs entirely and pushes the shared latent representation towards the assumed prior distribution $\alpha$ rather than the more natural barycenter distribution.
	Finally, for $k>2$, AlignFlow would require $\binom{k}{2}$ adversarial terms where each term adds significant complexity to training the model.
	Continuous normalizing flows based on solving ordinary differential equations (ODE) have added transportation costs to regularize the mappings and reduce the computational burden \citep{onken2020ot,finlay2020rnode}.
	However, these require using expensive ODE solvers whereas we focus on efficient explicitly invertible functions.

	\section{ITERATIVE ALGORITHM FOR SYMMETRIC MONGE MAP PROBLEM}

	To summarize the developments so far, we have cast the problem of learning a unique parsimonious shared latent representation given multiple distributions to that of solving a symmetric Monge map, which is equivalent to the OT barycenter problem.
	We now turn to the computational question of how to tractably estimate both the barycenter and the associated OT Monge maps from the respective class distributions to the barycenter.
	As mentioned in the introduction, estimating OT maps and barycenter distributions are hard in general---especially when only samples from the distributions are known---, and prior deep methods require complex adversarial optimization to handle the pushforward constraint.
	Thus, our starting point is that for \emph{certain classes of distributions}, the symmetric Monge map problem can be solved easily using well-known results from OT theory---thereby avoiding the need for adversarial learning.
	We provide specific examples of these classes of distributions and their corresponding symmetric Monge map solutions in \autoref{sec:symmetric-monge-layers}.
	Unfortunately, these simple classes of distributions and corresponding Monge maps may not be a good fit to the given samples, especially those in modern ML settings.
	Therefore, we propose an iterative approach that leverages the solutions to these simple problems at each iteration to build up a complex map one layer at a time.

	At each iteration, we first estimate the respective class densities $\{\hat{\mu_j}\}_{j=1}^{k}$ (possibly with shared parameters) from some carefully chosen simple classes of distributions (e.g., independent, Gaussian, density tree).
	We then estimate the corresponding symmetric Monge maps, denoted by $\{t_j\}_{j=1}^{k}$, which can be computed via known OT results for these simple classes of distributions.
	Thus, at iteration $M$, the transformed data correspond to transforms of the original distributions by the \textbf{deep maps} $T_j = t_j^{(M)} \circ t_j^{(M-1)} \circ \cdots \circ t_j^{(1)}$.
	We summarize this meta-algorithm in \autoref{alg:iterative-algorithm} and provide three layer examples in \autoref{sec:symmetric-monge-layers}.

	\begin{algorithm}[!ht]
		\caption{Iterative Symmetric Monge Meta-Algorithm}
		\label{alg:iterative-algorithm}
		\begin{algorithmic}
			\REQUIRE Samples from the $k$ class distributions $\mathcal{\xmat} = (X_1, X_2, \cdots, X_k)$, weight vector $\bm{\weight}$, number of iterations/layers $M$
			\ENSURE Estimated invertible symmetric Monge maps $(\text{T}} \newcommand{\supp}{\text{supp}_1, \text{T}} \newcommand{\supp}{\text{supp}_2, \cdots, \text{T}} \newcommand{\supp}{\text{supp}_k)$
			\STATE $\text{T}} \newcommand{\supp}{\text{supp}_j^{(0)} \gets \Id,\quad \forall j$
			\FOR{$\ell = \{1, 2, \dots, M\}$}
			\STATE $\hat{\mu}_1, \hat{\mu}_2, \cdots, \hat{\mu}_k \gets \textnormal{EstimateClassDensities}(\mathcal{X})$
			\STATE $t_1, t_2, \cdots, t_k \gets \textnormal{SymmetricMonge}(\hat{\mu}_1, \cdots, \hat{\mu}_k)$
			\FOR{$j = \{1, 2, \dots, k\}$}
			\STATE $X_j \gets t_j(X_j)$
			\STATE $\text{T}} \newcommand{\supp}{\text{supp}_j^{(\ell)} \gets t_j \circ \text{T}} \newcommand{\supp}{\text{supp}_j^{(\ell-1)}$
			\ENDFOR
			\ENDFOR

			\textbf{return} $(\text{T}} \newcommand{\supp}{\text{supp}^{(M)}_1, \text{T}} \newcommand{\supp}{\text{supp}^{(M)}_2, \cdots, \text{T}} \newcommand{\supp}{\text{supp}^{(M)}_k)$
		\end{algorithmic}
	\end{algorithm}

	\subsection{Related Iterative Approaches}
	Iterative Gaussianization  is an iterative density estimation method, that learns invertible flow-based models \citep{Chen2000, Lin2000, Lyu2009, Laparra2011, Balle2016}.
	The key idea is to first learn a rotation matrix via ICA \citep{ica1984} or similar method to linearly transform the data, and then Gaussianize each marginal independently.
	\citet{inouye2018deep} extend this by iteratively building normalizing flows from more general ``shallow'' density estimation approaches.
	More rencently, Sliced Iterative Generator (SIG) \citep{dai2020sliced} iteratively estimate an orthogonal linear transformation via sliced Wasserstein distance and then compute optimal transport on each dimension independently.
	However, these prior iterative approaches are focused on density estimation (i.e., learning a generative model), and in particular, learn a mapping between a \emph{known} base distribution (e.g. Gaussian) and the unknown data distribution.
	Whereas, our approach is focused on estimating a mapping between multiple \emph{unknown} distributions to a shared latent representation.

	Iterative approaches for mapping between two \emph{unknown} distributions include Projection Pursuit Monge Map \citep{meng2019ppmm} that iteratively finds interesting directions to project the data onto, and estimates Monge maps for the 1D projected data.
	The caveat, however, is that it used fixed interestingness functions such as variance to find the projection directions.
	\citet{kuang2019sample} propose an alternative iterative method for learning optimal maps and barycenters where each iteration requires the solution of a simpler but \emph{joint} optimal transport problem---rather than solving 1D OT problems as in \citep{meng2019ppmm}.
	In practice, \citet{kuang2019sample} use a set of fixed interestingness functions to find the needed structure.
	\citet{essid2019adaptive} extend this iterative approach by using an adversarial objective to automatically learn these interestingness functions.
	Our approach is more general since it does not require fixed interestingness functions or adversarial training.
	%

	\section{SYMMETRIC MONGE LAYERS}
	\label{sec:symmetric-monge-layers}
	We will now develop three examples of the symmetric Monge map layers that can be used in our meta-algorithm:
	a) multivariate Gaussian layer,
	b) independent components layer,
	and c) decision tree layer.
	For each layer, we detail the two key inner-loop steps from \autoref{alg:iterative-algorithm}: 1) estimating the class distributions given samples and 2) solving the symmetric Monge map given these class distributions.

	\subsection{Multivariate Gaussian Layer}

	\mysubsubsection{Multivariate Gaussian Distribution Estimators}
	We merely estimate a Gaussian distribution for each class \emph{independently} .
	The choice of the multivariate Gaussian estimator depends on the sampling regimes. In the small sample regime, we could choose a sparse Gaussian estimator such as graphical Lasso \citep{Friedman2008}.
	If we suspect some data corruption, we could use robust Gaussian estimators such as \citep{cheng2019faster}.
	With enough clean samples, we can merely use the sample mean and sample covariance for each class distribution, which we use in our experiments.

	\mysubsubsection{Multivariate Gaussian Symmetric Monge Maps}
	Given multivariate Gaussian class distributions, we can leverage known OT results to directly solve the symmetric Monge map problem (proof in appendix).
	{
		\renewcommand{\mu}{m}
		\renewcommand{\Sigma}{C}
		\begin{theorem}[Optimal symmetric Monge map for Gaussian distributions]
			\label{thm:gaussian-symmetric-monge}
			Given $k$ Gaussian distributions defined by their means and covariances $\{(\mu_j, \Sigma_j)\}_{j=1}^k$ and a weight vector $\bm{\weight}$, the optimal symmetric Monge map for Gaussian distributions is:
			$    T^*_j(\bm{\x}) \equiv T^*_{j \to \textnormal{bary}}(\bm{\x}) = \mu_{\textnormal{bary}} + A_{j\to\textnormal{bary}}(\bm{\x} - \mu_j) \, ,$
			where $A_{j \to \textnormal{bary}} \triangleq \Sigma_{\textnormal{bary}}^{-{\textstyle \frac{1}{2}}}\big(\Sigma_{\textnormal{bary}}^{{\textstyle \frac{1}{2}}} \Sigma_j \Sigma_{\textnormal{bary}}^{{\textstyle \frac{1}{2}}}\big)^{{\textstyle \frac{1}{2}}} \Sigma_{\textnormal{bary}}^{-{\textstyle \frac{1}{2}}}$, $\mu_{\textnormal{bary}}=\sum_{j=1}^k w_j \mu_j$, and $\Sigma_{\textnormal{bary}}$ is the unique fixed point of the map $\Psi(\Sigma) \triangleq \sum_{j=1}^k \sum_{j=1}^k w_j (\Sigma^{{\textstyle \frac{1}{2}}} \Sigma_j \Sigma^{{\textstyle \frac{1}{2}}})^{{\textstyle \frac{1}{2}}}$.
		\end{theorem}
	}

	\subsection{Independent Components (Na\"ive) Barycenter (NB)}
	\label{sec:independent-components}

	\mysubsubsection{Independent Component Distribution Estimators}
	We propose to estimate independent class distributions with an optional shared orthogonal transformation---similar to the independent components in independent component analysis (ICA) \citep{ica1984}.
	More formally, we can define our independent component class distributions as:
	$    \forall j \in \{1,\cdots,k\}$, $\mu_j(\bm{\x}; Q, \theta_j) = \prod_{s=1}^{\nvar} \mu_{js}([Q\bm{\x}]_{s}; \theta_{j,s})$,
	where $Q$ is an orthogonal matrix shared between class distributions, $\theta_j$ are the parameters for the $j$-th class distribution, and $[Q\bm{\x}]_{s}$ is the $s$-th coordinate after transforming $\bm{\x}$ by $Q$.
	While in theory we could optimize all class distributions jointly, if $Q$ is known, then the problem decouples into independent subproblems that can be solved very efficiently.
	Thus, we propose to first estimate the orthogonal transform $Q$ (which is the more challenging part) and then estimate the independent class distributions assuming that $Q$ is fixed.

	Although we could generate a random orthogonal matrix for $Q$, most random directions---particularly in high-dimensional spaces---will not reveal useful structure, which will yield slow progress towards a shared representation.
	Thus, we seek a $Q$ that reveals interesting structure between class distributions.
	We adopt the non-parametric approach of SIG that finds the orthogonal directions (or slices) that maximize the sum of 1D empirical Wasserstein distances along each direction \citep{dai2020sliced}.
	However, the formulation in \citet{dai2020sliced} only considers $k=2$ distributions.
	Thus, we generalize the approach to $k\geq 2$ distributions by minimizing the empirical \emph{symmetric} Wasserstein distances inspired by our symmetric problem:
	\begin{align}
	\argmax_{Q: Q^T Q = I_m} \!\! \bigg(\sum_{j=1}^k \!\frac{w_k}{m} \!\sum_{\ell=1}^m \!\frac{1}{n} \!\sum_{i=1}^n \!| (X_jq_\ell)_{[i]} - Y_{[i], \ell}|^p \!\bigg)^{\!\!\frac{1}{p}} \!\!\!,
	\label{eqn:mswd-distance}
	\end{align}
	where $p\geq 1$, $X_j$ is the sample data matrix for the $j$-th class, $(X_jq_\ell)_{[i]}$ signify the samples from the $j$-th class distribution projected along the direction $q_\ell$ sorted in ascending order, $Y_{[i], \ell} \triangleq \sum_{j=1}^k w_k (X_jq_\ell)_{[i]}$ is the empirical barycenter along direction $q_{\ell}$, $m \leq \nvar$ is the number of directions, and $I_m \in \ensuremath{\mathbb{R}}^{m \times m}$ is the identity matrix.
	For simplicity of exposition, we will assume that $m=\nvar$ for the rest of \autoref{sec:symmetric-monge-layers} but give details for $m < \nvar$ in the appendix.
	Intuitively, this finds the directions that reveal the largest difference between class distributions along each 1D projection.
	Importantly, note that our optimization for $Q$ naturally ignores any subspaces that have the same distribution among all classes---even if this shared distribution is arbitrarily complex; this yields more parsimonious transformations and preserves any shared structure in the class distributions.
	We adopt the approach of SIG that optimizes $Q$ directly on the manifold of orthonormal matrices (also called a Stiefel manifold) using projected gradient descent with backtracking line search (details in \citet{dai2020sliced}).

	Given the shared orthogonal matrix $Q$, we can now independently estimate the class distributions after projection, i.e., the class distributions of $\bm{\z}_j \equiv Q\bm{\x}_j$ for all $j\in\{1,\cdots,k\}$.
	Because we assume the class distributions are independent, we are free to use any univariate density estimators for each dimension independently including non-parametric estimators such as a Gaussian mixture or kernel density estimator.
	For high flexibility yet low computational cost, we choose to use a regularized histogram-based density estimator in our experiments (details in the appendix).
	We expect that any standard non-parametric estimator would work well.

	\mysubsubsection{Independent Component Symmetric Monge Maps}
	Given that the estimated class distributions are independent, we can leverage 1D OT results to efficiently compute the symmetric Monge map of each feature independently.
	\begin{theorem}[Optimal 1D Symmetric Monge Maps]
		\label{thm:univariate-symmetric-monge}
		The optimal univariate symmetric Monge maps are:
		$    T^*_j = F_{\textnormal{bary}}^{-1} \circ F_j \, ,$
		where $F_j$ is the CDF function of the $\mu_j$ distribution and $F_{\textnormal{bary}}^{-1}$ is the inverse CDF of the barycenter distribution, which is known to have the following form $F_{\textnormal{bary}}^{-1}(u) = \sum_j w_j F_j^{-1}(u)$.
	\end{theorem}
	To account for the shared orthogonal transformation $Q$, we can merely pre and post process the input by $Q$ and $Q^T$ respectively.
	Thus, the joint symmetric Monge map for the $j$-th class can be written as:
	$    T_j^* = Q^T \circ [T_{j,1}^*, \cdots, T_{j,\nvar}^*]^T \circ Q \,,$
	where $[T_{j,1}^*, \cdots, T_{j,\nvar}^*]^T$ denotes the component-wise transformation of each coordinate and $T_{j,s}$ is the $j$-th symmetric Monge map for dimension $s$.

	\subsection{Decision Tree Layer}
	To highlight that our framework can be used without gradient-based methods, we develop a decision tree layer.

	\mysubsubsection{Tree Distribution Estimators}
	Tree distributions are essentially multivariate generalizations of histograms with unequal bin sizes and splits based on the tree structure \citep{Ram2011,inouye2018deep}.
	Each leaf node has a constant density value and the overall tree density is thus piecewise constant.
	For this layer, we will assume the data has support on the unit hypercube.\footnote{
		If the data is not supported on the unit hypercube, the support can be modified by applying a simple preprocessor and postprocessor; for example, if the data lies in the real space, the CDF and inverse CDF of a standard normal distribution can be applied as a preprocessor and postprocessor respectively, which is what we do in our 2D experiments.
	}

	Similar to our independent components layer, our class distributions will share the same decision tree structure, but could have different density values for the leaves, i.e., the $j$-th class distribution has the form:
	$    p_j(\bm{\x}) = \{c_{j,\ell}, \text{ if } \bm{\x} \in \mathcal{N}_\ell\}$,
	where $\mathcal{N}_\ell$ is the $\ell$-th leaf node (i.e., the shared leaf node structure) and $c_{j,\ell}$ is the constant density value for that leaf node.
	Again, we choose to first estimate the shared tree structure $\mathcal{S}_{\textnormal{tree}}\triangleq \{\mathcal{N}_\ell\}_{\ell=1}^{L}$, and then estimate the density assuming this structure is fixed.
	To reveal the differences between class distributions, we use a decision tree classifier to find the shared tree structure.
	Note that the classifier will ignore any subspaces where the class distributions are the same (even if they are arbitrarily complex) similar to the independent layer; this will encourage more parsimonious maps and a preservation of shared structure if possible.
	Given the shared tree structure $\mathcal{S}_{\textnormal{tree}}$, the density values at each leaf node $c_{j,\ell}$ can be estimated via the smoothed leaf node frequencies as in \citet{inouye2018deep}.

	\mysubsubsection{Recursive Symmetric Monge Maps for Tree Densities}
	%
	Given the shared tree structure, we propose to estimate the symmetric Monge map by solving a sequence of 1D symmetric Monge map problems for all internal nodes in a bottom-up fashion.
	First, let us consider the base case of a single split.
	In this case, each class distribution can be viewed as a 2-bin histogram along the split dimension (while all other dimensions are assumed to be independent and uniform).
	Thus, we only need to solve the symmetric Monge problem for 1D histograms, which has a very efficient and closed-form solution by \autoref{thm:univariate-symmetric-monge} and the theorem below for barycenters of histograms (proof in appendix).
	To extend to full trees, we recursively solve these 2-bin histogram problems in a bottom-up fashion where the weight vector $\bm{\weight}$ is proportional to the class-specific probability for the node.
	\begin{theorem}[Barycenter of Histograms is a Histogram]
		\label{thm:barycenter-histogram}
		Given an ordered set of 1D histograms $(\mu_{\textnormal{hist},1}, \cdots, \mu_{\textnormal{hist},k})$ with a corresponding weight vector $\bm{\weight}$, the barycenter is another a histogram with at most the number of bin edges equal to the total number of bin edges across all histograms.
	\end{theorem}

	\section{EXPERIMENTS}
	We explore our iterative barycenter method both qualitatively and quantitatively using both 2D simulated data and ``permuted'' MNIST \citep{lecun-mnisthandwrittendigit-2010}---permuted means that our methods do not leverage the image structure of MNIST but merely treat each image as 784-dimensional vector.\footnote{AlignFlow and the CNN classifier we use for evaluation do use the image structure but our iterative methods do not.}  Additional experiments, implementation details and results can be found in the appendix, including experiments on FashionMNIST \citep{xiao2017/online}.

	\mysubsubsection{Metrics}
	We use pushforward constraint metrics to compare fairly with prior methods, but show transportation cost to highlight that our algorithm attempts to minimize it while prior methods ignore it.
	We estimate the transportation cost by an average over the test set, i.e., $\text{TC}=\sum_{j=1}^k \frac{w_j}{n_j} \sum_{x\in X_j} \|x- \hat{T}_j(x)\|^2$, where $X_j$ is the test dataset for the $j$-th class and lower is better.
	While transportation cost is simple to compute, it is challenging to determine how well the pushforward constraint (i.e., $T_{j\sharp}\mu_j = T_{j'\sharp}\mu_{j'},\quad \forall j\neq j'$) is satisfied---a critical piece for the symmetric Monge map problem.
	The key issue is that we only have samples from the class distributions $\mu_j$.
	We first note that the pushforward constraint can equivalently be written as $\mu_j =(T_{j}^{-1} \circ T_{j'})_{\sharp} \mu_{j'},\quad \forall j\neq j'$.
	Thus, for every class distribution $\mu_{j}$, we can sample $k-1$ ``fake'' distributions using our invertible transformations $\hat{\mu}_{j'\to j} = (\hat{T}_{j}^{-1} \circ \hat{T}_{j'})_{\sharp} \mu_{j'}$.
	For the 2D data, we can average the Wasserstein distance on test samples between the real and the fake distributions, i.e., $\text{WD} = \frac{1}{k^2 - k} \sum_{j\neq j'} \hat{W}(\mu_j, \hat{\mu}_{j'\to j})$, where $\hat{W}$ is the Wasserstein distance estimated using samples via the Sinkhorn algorithm with $\epsilon=0.1$ with maximum iterations set to 100.
	For higher dimensional data (e.g., MNIST), the Wasserstein distance between samples is a poor estimator of the true Wasserstein distance \citep{genevay2019sample}.
	Thus, we train a CNN classifier and compute the cross entropy loss between the real distribution $\mu_j$ and a uniform mixture over fake distributions $\bar{\mu}_{j'\to j} = \frac{1}{k-1}\sum_{j'} \hat{\mu}_{j'\to j}$ as proxy for measuring the distance between distributions, i.e., $\text{CEL} = \frac{1}{k} \sum_j \text{CELoss}(\mu_j, \bar{\mu}_{j'\to j}) = \frac{1}{k} \sum_j \min_f \E_{\mu_j}[-\log f(x)] + \E_{\bar{\mu}_{j'\to j}}[-\log (1-f(\bm{\x}))]$.\footnote{
		While in theory we could compute the cross entropy loss for every pair of class distributions, we do not consider that here because of the computational cost---each cross entropy loss term requires training a CNN classifier.}
	Higher CEL means it is harder to distinguish the real from the fake so higher is better.
	More details in the appendix.
	We compute the mean and standard deviation over 5 runs of each method.
	We also track approximate wall-clock training time for MNIST (all models trained on a CPU except AlignFlow which is trained on a GPU).

	\mysubsubsection{Methods}
	Given that there are no prior iterative approaches that map to a shared representation, we suggest a simple approach --- if a normalizing flow is estimated for each class distribution \emph{independently}, then the flows for each class can be seen to map to the same shared prior distribution.
	For this approach, the mappings and representations are not unique even in theory, and the shared distribution is not learned but merely assumed to be the prior distribution.
	However, this is the closest comparison that seems reasonable.
	Using this simple approach, we can compare our iterative algorithm to the flow-based iterative density estimation methods deep density destructors (DD) \citep{inouye2018deep} and sliced iterative generator (SIG) \citep{dai2020sliced}.
	For MNIST, as a non-iterative baseline, we compare to the invertible AlignFlow \citep{grover2020alignflow}, which explicitly maps both distributions to an assumed prior distribution.
	For our methods, we consider a single Gaussian barycenter (GB), a single independent (na\"ive) barycenter (NB), multi-layer tree barycenter (Tree), multi-layer random rotation followed by NB (Rand-NB), multi-layer orthogonal transform found using max sliced Wasserstein distance (\autoref{sec:independent-components}) followed by NB (mSWD-NB).
	Number of layers and other parameters are in the appendix.

	\begin{figure*}[!t]
		\newcommand{0.12\textwidth}{0.12\textwidth}
		\centering
		\begin{subfigure}[b]{0.16\textwidth}
			\includegraphics[width=\textwidth]{figures/two_toy/moon_ori.png}
			\caption{Original Data}
		\end{subfigure}
		\begin{subfigure}[t]{0.12\textwidth}
			\centering
			\includegraphics[width=\textwidth]{figures/two_toy/moon_gb.png}
			\vspace{-2em}
			\caption{GB}
		\end{subfigure}
		\begin{subfigure}[t]{0.12\textwidth}
			\centering
			\includegraphics[width=\textwidth]{figures/two_toy/moon_swdnb.png}
			\vspace{-2em}
			\caption{mSWD-NB}
		\end{subfigure}
		\begin{subfigure}[t]{0.12\textwidth}
			\centering
			\includegraphics[width=\textwidth]{figures/two_toy/moon_tree.png}
			\vspace{-2em}
			\caption{Tree}
		\end{subfigure}
		\begin{subfigure}[t]{0.12\textwidth}
			\centering
			\includegraphics[width=\textwidth]{figures/two_toy/moon_dd.png}
			\vspace{-2em}
			\caption{DD}
		\end{subfigure}
		\begin{subfigure}[t]{0.3\textwidth}
			\includegraphics[width=\textwidth]{figures/two_toy/wd_test_0.2.png}
			\vspace{-2em}
			\caption{WD over 30 Layers}
			\label{fig:convergence}
		\end{subfigure}
		\vspace{-0.5em}
		\caption{The shared representations (top row) for each method show that our iterative methods (DD, Tree) find low transportation cost shared latent spaces whereas DD ignores transportation cost and merely projects both distributions to the uniform distribution.
			The bottom row shows test samples that were flipped to the other class distribution (ideally these ``fake'' samples would look like the original data).
		}
		\label{fig:2d-k2-qualitative}
	\end{figure*}

	\begin{figure*}[!t]
		\centering
		\begin{subfigure}[t]{0.24\textwidth}
			\centering
			\includegraphics[width=\textwidth]{figures/two_mnist/gb0.png}
			\caption{GB}
		\end{subfigure}
		\begin{subfigure}[t]{0.24\textwidth}
			\centering
			\includegraphics[width=\textwidth]{figures/two_mnist/swdnb0.png}
			\caption{mSWD-NB}
		\end{subfigure}
		\begin{subfigure}[t]{0.24\textwidth}
			\centering
			\includegraphics[width=\textwidth]{figures/two_mnist/dd0.png}
			\caption{DD}
		\end{subfigure}
		\begin{subfigure}[t]{0.24\textwidth}
			\centering
			\includegraphics[width=\textwidth]{figures/two_mnist/alignflow0.png}
			\caption{AlignFlow}
		\end{subfigure}
		\vspace{-1em}
		\caption{These examples demonstrate that our methods find a more natural shared latent representation that preserves structural similarities (e.g., black pixels) between the two digits while DD and AlignFlow do not. The rows from top to bottom are the original MNIST digits, their latent representation, and their projection to the space of the other digit (i.e., flipped).}
		\label{fig:mnist-k2-samples}
		\vspace{-1em}
	\end{figure*}
	\mysubsubsection{2D experiments}
	The qualitative results in \autoref{fig:2d-k2-qualitative} illustrate that our methods (Tree, mSWD-NB) find shared latent spaces closer to the barycenter distribution (since the transportation cost is low) whereas density destructors (DD) ignores transportation costs and projects both distributions to the uniform distribution.
	The results for the 2D datasets with $k=2$ in \autoref{tab:2d-k2-barycenter-flow} demonstrate that our iterative barycenter flows perform comparable or better than the baseline methods (DD, SIG) in terms of the pushforward constraint, which is measured by the empirical Wasserstein-2 distance on test data (WD), while having significantly lower transportation cost (TC) on test data.
	Additional experiments and results for $k>2$ in appendix.
	\begin{table}[!ht]
		\vspace{-1em}
		\centering
		\caption{
			Transportation cost (TC), sample-based Wasserstein distance (WD, lower is better) for 2D data.
			More 2D datasets in appendix.
		}
		\label{tab:2d-k2-barycenter-flow}
		\includegraphics[width=0.4\linewidth]{tables/2d-k2-moons-barycenter-flow-results.png}
		\vspace{-1em}
	\end{table}

	\mysubsubsection{``Permuted'' MNIST}
	Qualitative samples from the latent space and after flipping between the two digits (\autoref{fig:mnist-k2-samples}) highlight that our barycenter methods retain shared latent structure such as the black pixels, whereas the generative baselines (DD, AlignFlow) try to move the shared latent distribution to the assumed prior (uniform or Gaussian, respectively) so that shared structure is also removed.
	To emphasize that our method can scale to more class distributions, we present qualitative examples of transforming between every digit and every other digit (i.e., $k=10$) for MNIST in \autoref{fig:mnist-10-samples} (quantitative results in appendix).
	In \autoref{fig:mnist-10-samples}, real samples are along the diagonal and each column shows the mapping to another digit distribution (e.g., the bottom left ``9'' image is the mapping of the real sample ``0'' on the top left).
	Notice that even for this multiclass case, almost all transformed digits are recognizable.
	This multi-class situation (i.e. $k>2$) is much more difficult for AlignFlow (which did not implement $k>2$) and would na\"ively require $k^2 - k$ pairwise adversarial loss terms.
	We emphasize that our extension to $k>2$ is straightforward because we only need to solve simple symmetric Monge map problems, which scale linearly with $k$.
	The results for $k=2$ in \autoref{tab:MNIST-k2} show that our methods perform well in terms of the pushforward constraint (measured by the CNN cross-entropy loss where higher is better) while having significantly lower transportation costs (TC) than the iterative baseline (DD) and end-to-end baseline AlignFlow.
	Also, the computational cost is much lower for the iterative methods ($<1.5$ hour on CPU), whereas AlignFlow trained for 50 epochs on a GPU took approximately 11 hours (thus, we only estimate one model and cannot compute standard deviations for AlignFlow).
	While the tree approach performed comparably in low dimensions, the tree approach performs poorly in this high-dimensional setting.
	\begin{table}[!t]
		\centering
		\vspace{-1em}
		\caption{
			Transportation cost (TC), cross entropy loss (CEL, higher is better) and time for MNIST.
		}
		\label{tab:MNIST-k2}
		\includegraphics[width=0.4\linewidth]{tables/MNIST-k2-barycenter-flow-results.png}
		\vspace{-1em}
	\end{table}
	\begin{figure}[!ht]
		\centering
		\includegraphics[width=0.75\linewidth]{figures/ten_mnist/swdnb0.png}
		\vspace{-0.5em}
		\caption{Multi-domain ($k=10$) results for MNIST with mSWD-NB.
			The diagonal shows the real samples and each column shows the mappings of the real sample to the distribution of the other digits.}
		\label{fig:mnist-10-samples}
		\vspace{-1.5em}
	\end{figure}

	\section{DISCUSSION AND CONCLUSION}

	We seek to estimate the maps between two or more distributions and a unique and parsimonious shared latent representation.
	We formulate the symmetric Monge map problem and prove that it maps to the natural barycenter distribution---a result that may be of independent interest.
	Unlike normalizing flows or standard generative approaches, our formulation does not require a pre-fixed latent distribution but can be applied to any two \emph{unknown} distributions.
	However, this requires solving a very challenging OT problem both theoretically and computationally.
	Therefore, we leverage OT theory and insights from modern density estimation to construct a practical iterative estimation algorithm.
	Specifically, each iteration of our algorithm adds a layer to our model by first estimating the class distributions and then estimating the symmetric Monge map for this problem using known OT results.
	Unlike prior methods, our approach enables efficient \emph{multi-class} distribution mappings to a shared representation (e.g., see \autoref{fig:mnist-10-samples}).
	More generally, our iterative approach can also provide a complementary perspective for solving challenging deep learning problems because it is computationally efficient and does not require end-to-end learning, gradient-based methods, or adversarial optimization.
	In our experiments, we see that the shared latent space indeed seems qualitatively better than alternatives because it preserves some shared structure (\autoref{fig:mnist-k2-samples}).

	We believe our work is a first step that can open up novel perspectives, algorithms, and theory.
	For example, we expect that the sample complexity of our symmetric Monge problem will be lower than normalizing flow approaches such as AlignFlow or density destructors.
	To build intuition, consider that images lie on lower dimensional curves.
	Normalizing flows would need to transform these curves to a Gaussian or uniform distribution over the entire ambient image space.
	However, the symmetric Monge problem would only have to transform both distributions to some shared barycenter curves.
	If the images share structure, these curves will be much closer to each other than a Gaussian distribution.
	We leave development of this theory to future work.

	Despite many advantages of our approach, however, there are also many open challenges.
	For example, can we guarantee that each iteration of our algorithm reduces some statistical distance between the distributions similar to what has been proved for iterative Gaussianization?
	While we initially attempted to prove this, we believe theoretical guarantees of our iterative algorithm to be much more challenging than iterative Gaussianization because the class distributions are assumed to be \emph{unknown} except through samples (while Gaussianization assumes one distribution is a known prior distribution).
	Given the generality of the class distributions, we expect that while counterexamples likely exist (thereby excluding simple proofs), these worst cases are uncommon in practice given our experimental results.
	We expect that modern smoothed analysis could enable proofs of convergence beyond simple worst-case analysis, but this will require significant theoretical development so we leave this to future work.
	Another theoretical challenge is proving approximation bounds on how well the iterative algorithm approximates the true barycenter.
	While our experiments demonstrate that the approximation error is reasonable in practice, providing bounds for our iterative algorithm even under certain regularity conditions of the unknown class distributions would be an excellent future direction.
	Thus, we believe that our formulation as a symmetric Monge problem and our iterative algorithm open up new theoretical and fundamental questions.

	\section*{Acknowledgments}
	D.I., Z.Z., and Z.G. acknowledge support from the Army Research Lab through contract number W911NF-2020-221. P. R. acknowledges the support of NSF via IIS-1909816.
	\clearpage
	{
		\bibliographystyle{plainnat}


\section{The tight-binding Hamiltonian of the model system}
\label{S-1}

\begin{figure}
\centering
\includegraphics[width=0.7\textwidth]{Figs1.pdf}
\caption{
  (a)  The spatial distribution of the confining potential $V^{}_{\rm DQD-in}(x)$ that determines the intra-cell tunneling amplitude within the double quantum-dot chain in the main text (see Fig.~1). (b) Confining potential $V^{}_{\rm DQD-ex}(x)$ that determines the inter-cell tunneling amplitude. (c) The intra/inter-cell hopping amplitude of an InAs double quantum dot chain, $t^{}_{\rm in/ex,0}$, as a function  of the potential difference $\Delta V=V^{}_{1}-V^{}_{2}$, with the  periodic length $2L=120$~nm and  $d=10$~nm. Here, $t^{}_{\rm in,0}$ ($t^{}_{\rm ex,0}$) is denoted by the solid (dashed) line and, meanwhile, the variations of the  chemical potential $V_{\rm c}$ are also displayed by the dash-doted curve.  (d)   The intra- and inter-cell spin-flip tunneling amplitudes, $t^{\prime}_{\rm in}$ and  $t^{\prime}_{\rm ex}$, of the double-quantum-dot chain as functions of the tilting angle $\theta$, with the potential difference $\Delta V^{}_{}=-10$~mev,  the  magnetic field strength $B=1.0$~T, and the spin-orbit coupling scaled by $x^{}_{\rm so}\equiv \hbar/(m_{e}\alpha)=180$~nm.}
\label{Figs1}
\end{figure}

 Here is outlined the  derivation of the  nearest-neighbor tight-binding (TB) Hamiltonian of our model.
 The derivation is based on the original Hamiltonian, Eq.~(5) in the main text. The periodic gate potential there, $U(x)$,  is viewed as composed of two alternating quantum dots, QD1 and QD2, as depicted in Fig.~\ref{Figs1}(a) and (b). Each quantum dot obeys the Hamiltonian
 \begin{align}
 H^{}_{\rm 1/2}=\frac{p^{2}}{2m_{e}}+ V^{}_{\rm qd,1/2}(x)+\alpha p\sigma_{y}+\frac{\Delta_{\rm z}}{2}\sigma^{}_{n},
 \label{H12}
 \end{align}
 where  $\alpha$ is  the strength of the Rashba spin-orbit coupling (SOC), $\Delta_{\rm z}$ is the Zeeman splitting, $\sigma^{}_{n}= \hat{\boldsymbol{\sigma}}\cdot \hat{\mathbf{n}}$ with  $\hat{\mathbf{n}}=\{0, \cos\theta, \sin\theta\}$ being the direction vector of the external magnetic field  and  $\hat{\boldsymbol{\sigma}}$ denotes  the Pauli matrices,  $\hat{\boldsymbol{\sigma}}=\{\sigma^{}_{x},\sigma^{}_{y},\sigma^{}_{z}\}$.  The  confining  potentials that define the quantum dots are
 \begin{align}
 V^{}_{\rm qd,1}(x)=\begin{cases}
 V_{1} ~~~~~~~~~~~~~x\geq -d\\
 0~~~~~~~ -L+d\leq x<-d,\\
 V_{2}~~~~~~~~~~~x<-L+d
 \end{cases}~~~~   V^{}_{\rm qd,2}(x)=\begin{cases}
 V_{1} ~~~~~~~~~~~~~x< d\\
 0~~~~~~~ d\leq x<L-d.\\
 V_{2}~~~~~~~~x\geq L-d
 \end{cases}
 \label{lcp}
 \end{align}

 We next solve for the energy eigenstates of the two dots using Eq. (\ref{H12}). Moreover,  in order to preserve the locality of the wave functions, we assume  the electron energy eigenvalue   $E<V^{}_{1},V^{}_{2}$. It
is expedient to consider the Schr\"{o}dinger equation pertaining to all three regions in Eq. (\ref{lcp}) by introducing
 \begin{align}
 \widetilde{H}^{}_{\zeta=0,1,2}\psi^{}_{\zeta }(x)=E\psi_{\zeta}(x).
 \label{WQD}
 \end{align}
 Here,   $\psi^{}_{\zeta}$ is the corresponding eigenfunction on each portion  and
\begin{align}
\widetilde{H}^{}_{\zeta}=\frac{p^{2}}{2m_{e}}+\alpha p\sigma^{}_{y}+\frac{\Delta^{}_{\rm z}}{2}\sigma^{}_{n}+V^{}_{\zeta},
\label{LQD}
 \end{align}
  where
 $V_{0}=0$, and $V_{1,2}$ are given in Eqs. (\ref{lcp}).
The eigenfunctions of $\widetilde{H}^{}_{\zeta}$ are spinors,
 \begin{align}
 \psi^{}_{\zeta}(x)=\exp\left(i\gamma^{}_{\zeta}x\right)\begin{pmatrix}
 w_{\zeta}\\
 d_{\zeta}
 \end{pmatrix}.
 \label{BEFs}
 \end{align}
For a fixed value of $E$ there exist four independent solutions to Eq.~(\ref{WQD}), with the corresponding wave vectors  $\gamma^{}_{\zeta, \lambda=1,2,3,4}$  determined by the  quartic equation
  \begin{align}
  \left(\frac{\hbar^{2}{\gamma^{2}_{\zeta}}}{2m_{e}}+V^{}_{\zeta}-E\right)^{2}-\hbar^{2}\gamma^{2}_{\zeta}\alpha^{2}-\Delta_{\rm z}\hbar \gamma^{}_{\zeta}\alpha\cos(\theta)-\frac{\Delta^{2}_{\rm z}}{4}=0.
  \label{wk}
  \end{align}
  The respective spinor components are
 \begin{align}
w^{}_{\zeta,\lambda}= i\hbar\alpha\gamma^{}_{\zeta, \lambda}+i\frac{\Delta^{}_{\rm z}}{2}\cos(\theta),~~~~~~d^{}_{\zeta,\lambda}=\frac{\hbar^{2}\gamma^{2}_{\zeta, \lambda}}{2m_{e}}+\frac{\Delta_{\rm z}}{2}\sin(\theta)+V^{}_{\zeta}-E.
  \label{wkc}
 \end{align}
 In particular, when $\theta=\pi/2$, the expressions for the wave vectors can be derived analytically
 \begin{align}
 \gamma^{}_{\zeta,1/2}=&\pm\sqrt{\frac{2m_{e}}{\hbar^{2}}}\sqrt{E-V_{\chi}+m_{e}\alpha^{2}+\sqrt{m^{2}_{e}\alpha^{4}+2m_{e}\alpha^{2}(E-V_{\chi})+\Delta^{2}_{\rm z}/4}} \nonumber\\
 \gamma^{}_{\zeta,3/4}=&\pm\sqrt{\frac{2m_{e}}{\hbar^{2}}}\sqrt{E-V_{\chi}+m_{e}\alpha^{2}-\sqrt{m^{2}_{e}\alpha^{4}+2m_{e}\alpha^{2}(E-V_{\chi})+\Delta^{2}_{\rm z}/4}}.

 \end{align}
The convergence of the wave function at  $|x|\rightarrow \infty$ requires the wave vectors to have ${\rm Im}[\gamma^{}_{1}]>0$  for dot 1 and ${\rm Im}[\gamma^{}_{2}]>0$ for dot 2.  Meanwhile, the convergence condition when $x\rightarrow -\infty$ demands  ${\rm Im}[\gamma^{}_{2}]<0$  for dot 1 and ${\rm Im}[\gamma^{}_{1}]<0$ for dot 2.
Then, the energy eigenfunctions of QD1 and QD2 can be written as
\begin{align}
\Psi_{1}(x)=&\sum^{4}_{\lambda=1}c^{}_{0,\lambda}[1-\Theta(x+d)]\Theta(x+L-d)\psi^{}_{0,\lambda}(x)+\Theta(x+d)[c^{}_{1,1}\psi^{}_{1,1}(x)+c^{}_{1,4}\psi^{}_{1,4}(x)]\nonumber\\
&~~~~~~~+[1-\Theta(x+L-d)][c^{}_{2,2}\psi^{}_{2,2}(x)+c^{}_{2,3}\psi^{}_{2,3}(x)]\\
\Psi_{2}(x)=&\sum^{4}_{\lambda=1}v^{}_{0,\lambda}[1-\Theta(x-L+d)]\Theta(x-d)\psi^{}_{0,\lambda}(x)+[1-\Theta(x-d)][v^{}_{1,2}\psi^{}_{1,2}(x)+v^{}_{1,3}\psi^{}_{1,3}(x)]\nonumber\\
&~~~~~~~+\Theta(x-L+d)[v^{}_{2,1}\psi^{}_{2,1}(x)+v^{}_{2,4}\psi^{}_{2,4}(x)],
\end{align}
where $\Theta(x)$ is the Heaviside step function  $\Theta(x)=\begin{cases}
0  ~~x\leq0\\
1 ~~x>0\end{cases}$, the coefficients $c^{}_{\zeta,\lambda}$ and  $v^{}_{\zeta,\lambda}$  are determined by the boundary conditions, i.e., from the continuity  equations at $x^{}_{1/3}=\mp d$ and $x^{}_{2/4}=\mp L\pm d$.
 Using the explicit form for $V^{}_{qd,1}(x)$ in Eq.~(\ref{lcp}), the boundary conditions for the localized wave function on QD1 are \cite{com}
\begin{align}
\lim_{\epsilon \rightarrow 0 }[\Psi_{1}(x_{1}+\epsilon)-\Psi_{1}(x_{1}-\epsilon)]=0 ,~~~~\lim_{\epsilon\rightarrow 0} [\Psi_{1}(x_{2}+\epsilon)-\Psi_{1}(x_{2}-\epsilon)]=0 , \nonumber\\
\lim_{\epsilon \rightarrow 0 }[\Psi'_{1}(x_{1}+\epsilon)-\Psi'_{1}(x_{1}-\epsilon)]=0 ,~~~~\lim_{\epsilon\rightarrow 0} [\Psi'_{1}(x_{2}+\epsilon)-\Psi'_{1}(x_{2}-\epsilon)]=0.
\label{CD-1}
\end{align}
By regarding $\hat{z}_{}^{T}= \{c^{}_{0,1},c^{}_{0,2},c^{}_{0,3},c^{}_{0,4},c^{}_{1,1},c^{}_{1,4},c^{}_{2,2},c^{}_{2,3}\}$ as the variable vector, Eq. (\ref{CD-1}) can be written in a matrix form
\begin{align}
\mathbf{M} \cdot \hat{z}=0,
\label{MAE}
\end{align}
where
\begin{footnotesize}
\begin{align}
\mathbf{M} =\begin{pmatrix}
\psi^{}_{0,1,1}(x_{1})&\psi^{}_{0,2,1}(x_{1})&\psi^{}_{0,3,1}(x_{1})&\psi^{}_{0,4,1}(x_{1})&-\psi_{1,1,1}(x_{1})&-\psi_{1,4,1}(x_{1})&0&0\\
\psi^{}_{0,1,2}(x_{1})&\psi^{}_{0,2,2}(x_{1})&\psi^{}_{0,3,2}(x_{1})&\psi^{}_{0,4,2}(x_{1})&-\psi_{1,1,2}(x_{1})&-\psi_{1,4,2}(x_{1})&0&0\\
\psi^{}_{0,1,1}(x_{2})&\psi^{}_{0,2,1}(x_{2})&\psi^{}_{0,3,1}(x_{2})&\psi^{}_{0,4,1}(x_{2})&0&0&-\psi_{2,2,1}(x_{2})&-\psi_{2,3,1}(x_{2})\\
\psi^{}_{0,1,2}(x_{2})&\psi^{}_{0,2,2}(x_{2})&\psi^{}_{0,3,2}(x_{2})&\psi^{}_{0,4,2}(x_{2})&0&0&-\psi_{2,2,2}(x_{2})&-\psi_{2,3,2}(x_{2})\\
\gamma^{}_{0,1}\psi^{}_{0,1,1}(x_{1})&\gamma^{}_{0,2}\psi^{}_{0,2,1}(x_{1})&\gamma^{}_{0,3}\psi^{}_{0,3,1}(x_{1})&\gamma^{}_{0,4}\psi^{}_{0,4,1}(x_{1})&-\gamma^{}_{1,1}\psi_{1,1,1}(x_{1})&-\gamma^{}_{1,4}\psi_{1,4,1}(x_{1})&0&0\\
\gamma^{}_{0,1}\psi^{}_{0,1,2}(x_{1})&\gamma^{}_{0,2}\psi^{}_{0,2,2}(x_{1})&\gamma^{}_{0,3}\psi^{}_{0,3,2}(x_{1})&\gamma^{}_{0,4}\psi^{}_{0,4,2}(x_{1})&-\gamma^{}_{1,1}\psi_{1,1,2}(x_{1})&-\gamma^{}_{1,4}\psi_{1,4,2}(x_{1})&0&0\\
\gamma^{}_{0,1}\psi^{}_{0,1,1}(x_{2})&\gamma^{}_{0,2}\psi^{}_{0,2,1}(x_{2})&\gamma^{}_{0,3}\psi^{}_{0,3,1}(x_{2})&\gamma^{}_{0,4}\psi^{}_{0,4,1}(x_{2})&0&0&-\gamma^{}_{2,2}\psi_{2,2,1}(x_{2})&-\gamma^{}_{2,3}\psi_{2,3,1}(x_{2})\\
\gamma^{}_{0,1}\psi^{}_{0,1,2}(x_{2})&\gamma^{}_{0,2}\psi^{}_{0,2,2}(x_{2})&\gamma^{}_{0,3}\psi^{}_{0,3,2}(x_{2})&\gamma^{}_{0,4}\psi^{}_{0,4,2}(x_{2})&0&0&-\gamma^{}_{2,2}\psi_{2,2,2}(x_{2})&-\gamma^{}_{2,3}\psi_{2,3,2}(x_{2})
\end{pmatrix},
\label{MES}
\end{align}
\end{footnotesize}
with the two spin components $\psi^{}_{\zeta,\lambda,1}(x)=w^{}_{\zeta,\lambda}\exp(i\gamma^{}_{\zeta,\lambda} x)$ and $\psi^{}_{\zeta,\lambda,2}(x)= d^{}_{\zeta,\lambda}\exp(i\gamma^{}_{\zeta,\lambda} x)$. Obviously, the nontrivial solution to Eq.~(\ref{MES}) requires
\begin{align}
\det[ \mathbf{M}(E)]=0,
\label{DEM}
\end{align}
which yields an equation for the energy eigenvalues $E$, see Eqs.~(\ref{wk})  and (\ref{wkc}).
 Let $E_{ \Uparrow}$ and $E_{ \Downarrow}$  denote the two lowest energy solutions to Eq.~(\ref{DEM}), with $E_{ \Uparrow}\geq E_{ \Downarrow}$. The explicit expressions for the two corresponding spinors can also be completely determined by the boundary conditions   Eqs.~ (\ref{CD-1}) in conjunction with the normalization condition $\int \Psi^{\dagger}_{1,\tau=\Uparrow,\Downarrow}(x)\Psi^{}_{1,\tau}(x)dx=1$.
where $\tau=\Uparrow,~\Downarrow$ indicate the two quasi-spin eigenstates on the dot.

The two lowest-energy spin states of the other quantum dot, QD2,
i.e., $\Psi^{}_{2, \Uparrow}(x)$ and $\Psi^{}_{2,\Downarrow }(x)$,
 can be derived similarly to the discussion above. Due to the inversion symmetry  which implies  $V^{}_{\rm qd,1}(x)=V^{}_{\rm qd,2}(-x)$, the corresponding energy eigenvalues are identical to  those of QD1.

 Exploiting the localized wave functions on each quantum dot, we next derive the tunneling amplitudes of the double-quantum-dot (DQD) model, utilizing linear combinations of the atomic orbitals.
The intra-cell amplitude is expediently found by considering  the interdot tunneling in the double
 quantum-dot Hamiltonian [see Fig.~\ref{Figs1}(a)]
 \begin{align}
 H^{}_{\rm DQD-in}(x)=\frac{p^{2}}{2m_{e}}+\alpha p\sigma_{y}+\frac{\Delta^{}_{\rm z}}{2}\sigma_{n}+V^{}_{\rm DQD-in}(x)-V^{}_{\rm c},
 \label{ODQD1}
 \end{align}
with the double-well potential
\begin{align}
 V^{}_{\rm DQD-in}(x)=V^{}_{1}[1-\Theta(|x|-d)]+V^{}_{2}\Theta(|x|-L+d)
\end{align}
and the constant chemical potential $V^{}_{\rm c}=\frac{1}{2}(E^{}_{\Uparrow}+E^{}_{\Downarrow})$ to offset the on-site orbital energy (which is conventionally set at zero for the "standard" SSH model).
To this end, one first derives the orthonormal basis $\left\{\Phi^{}_{1,\Uparrow}(x),\Phi^{}_{1,\Downarrow}(x),\Phi^{}_{2,\Uparrow}(x),
\Phi^{}_{2,\Downarrow}(x)\right\}$ using the four localized states $\Psi^{}_{1/2,\Uparrow}(x)$ and  $\Psi^{}_{1/2,\Downarrow}(x)$. Those are then used to present the electron field operators as
$
\Psi_{\rm e }(x)= \sum^{}_{\sigma=\Uparrow,\Downarrow}[a^{}_{\sigma}\Phi^{}_{1,\sigma}(x) + b^{}_{\sigma}\Phi^{}_{2,\sigma}(x)]$,
 where $a^{}_{\sigma}$ ($a^{\dagger}_{\sigma}$) and $b^{}_{\sigma}$ ($b^{\dagger}_{\sigma}$) are the annihilation (creation) operators for the quasi-spin states on QD1 and QD2,
 respectively. The second quantized form of $H^{}_{\rm DQD-in}$ is then
 \begin{align}
 H ^{}_{\rm DQD-in}=& \int \Psi^{\ast}_{\rm e}(x) H^{}_{\rm DQD-in}(x) \Psi^{}_{\rm e }(x)dx \nonumber\\
=&\frac{\Delta^{}_{\rm z}}{2} \left(a^{\dagger}_{ }\sigma^{}_{z}a^{}_{ }+b^{\dagger}_{ }\sigma^{}_{z}b^{}_{ }\right)+\left(a^{\dagger}t^{}_{\rm in}b+{\rm H.c.}\right)\ ,
\label{DQD-in}
\end{align}
with  the spinors $a^{}_{}=\{a^{}_{\Uparrow},a^{}_{\Downarrow}\}$ and $b^{}_{}= \{b^{}_{\Uparrow},b^{}_{\Downarrow}\}$. In fact, the on-site energies of different quasi-spin states in each dot are modified by the spin-orbit interaction. However, as compared to the large spin-orbit length of the InAs nanowire, i.e., $x^{}_{\rm so}=180$~nm, the small characteristic length of the dots,  i.e., $L-2d=40$~nm, renders the SOC-induced modification negligible, and then it is a good approximation to write  the  quasi-spin up (down) on-site energy as the positive (negative) half Zeeman splitting energy as shown in Eq.~(\ref{DQD-in}). The quasi-spin tunneling amplitudes can be written as a matrix
\begin{align}
t^{}_{\rm in}=t^{}_{{\rm in}, 0}\exp\left[i\varphi^{}_{{\rm in} } (\hat{\mathbf{v}}^{}_{\rm in}\cdot\hat{\boldsymbol{\sigma}})\right]\ ,
\label{T-IN}
\end{align}
with
\begin{align}
[t^{}_{\rm in}]^{}_{\sigma \sigma'}=\int \Phi^{\dagger} _{1,\sigma}(x)H^{}_{\rm DQD-in}\Phi_{2,\sigma' }(x)dx\ .
\label{s19}
\end{align}

As seen, the effect of the  spin-orbit coupling appears as a (matrix) phase factor multiplying the intra-cell spin tunneling matrix.
 Here, $t^{}_{\rm in,0}$ is the real tunneling amplitude dominated by the potential difference $\Delta V=V_{1}-V_{2}$ as shown in Fig.~\ref{Figs1}(c), $\varphi^{}_{\rm in }$ is the spin rotation phase induced by the spin-orbit coupling, and $\hat{\mathbf{v}}^{}_{\rm in}=\{\sin(\vartheta^{}_{\rm in}),0,\cos(\vartheta^{}_{\rm in})\}$ is the unit vector determining the direction of the axis of rotation. Numerically, $\varphi^{}_{\rm in }$ is proportional to the ratio between the interdot distance and the spin-orbit length $L/x^{}_{\rm so}$ with $x^{}_{\rm so}= \hbar/(m^{}_{e}\alpha)$,  and that is consistent with
the analytical analysis in Ref.~\cite{Liu2018}

It is important to note that the coordinate system for defining $\hat{\mathbf{v}}^{}_{{\rm in}}$ spin state
 does not coincide with that defining the spin operators in Eq.~(\ref{ODQD1}).
The $x$-component of $\hat{\mathbf{v}}^{}_{{\rm in}}$ corresponds to the interdot spin-flip tunneling between the localized quasi-spin states, and vanishes in the absence of spin mixing.
Therefore, for zero magnetic field, i.e., when $\Delta^{}_{\rm z}=0$, or when $\theta=0$, the vector $\hat{\mathbf{v}}^{}_{\rm in}$ is purely along the $z$ direction since then $[H^{}_{\rm DQD-in},\sigma^{}_{y}]=0$. In other cases,
the specific value of $\vartheta^{}_{\rm in }$ depends on the degree of spin mixing on each localized quasi-spin state, which in turn is determined by the angle between the external magnetic field and the one induced by the spin-orbit  interaction. The amplitude of the spin-flip tunneling in Eq.~(\ref{T-IN}), $t^{\prime}_{\rm in }=t^{}_{\rm in, 0}\sin(\varphi^{}_{\rm in }) \sin(\vartheta^{}_{\rm in })$, as a function of the tilting angle $\theta$ for a fixed magnitude of the magnetic field, is portrayed in Fig.~\ref{Figs1}(d), showing that it vanishes for $\theta=0$.

The inter-cell hopping is addressed by considering a neighboring DQD [as shown in Fig.~\ref{Figs1}(b)] which obeys the Hamiltonian
\begin{align}
H^{}_{\rm DQD-ex}(x)=\frac{p^{2}}{2m_{e}}+\alpha p \sigma_{y}+\frac{\Delta_{\rm z}}{2}\sigma^{}_{n}+V^{}_{\rm DQD-ex}(x)-V^{}_{\rm c},
\end{align}
where the double-well potential takes the form
\begin{align}
V^{}_{\rm DQD-ex}(x)=V^{}_{2}[1-\Theta(|x|-d)]+V^{}_{1}\Theta(|x|-L+d).
\end{align}
Similar to the intra-cell spin tunneling discussed above, the spin transfer matrix for the inter-cell tunneling  takes the form
\begin{align}
t^{}_{\rm ex}= t^{}_{{\rm ex},0}\exp\left[i\varphi^{}_{{\rm ex} } (\hat{\mathbf{v}}^{}_{\rm ex}\cdot\hat{\boldsymbol{\sigma}})\right]\ ,
\label{T-EX}
\end{align}
with $t^{}_{{\rm ex},0}$ and $\varphi^{}_{{\rm ex}}$ corresponding to the real amplitude and the spin rotation phase of the interdot tunneling, and $\hat{\mathbf{v}}^{}_{{\rm ex}}=\{\sin(\vartheta^{}_{{\rm ex}}),0,\cos(\vartheta^{}_{{\rm ex}})\}$  describing the direction of the axis of rotation.
 Analogously, the corresponding spin-flip tunneling amplitude  $t^{\prime}_{\rm ex}=t^{}_{{\rm ex},0}\sin(\varphi^{}_{\rm ex})\sin(\vartheta^{}_{\rm ex})$ can be regulated by manipulating the tilting angle $\theta$, as shown in Fig.~\ref{Figs1}(d).
Therefore, by including the intra-cell spin tunneling, the second quantized Hamiltonian of the double-quantum-dot chain is
\begin{align}
H^{}_{\rm T}= \sum^{}_{n}\frac{\Delta^{}_{\rm z}}{2}\left(a^{\dagger}_{n} \sigma^{}_{z}a^{}_{n}+b^{\dagger}_{n}\sigma^{}_{z}b_{n} \right)+\sum^{}_{n}\left(a^{\dagger}_{n}t^{}_{\rm in }b^{}_{n}+a^{\dagger}_{n+1}t^{}_{\rm ex}b^{}_{n}+{\rm h.c.}\right)\ ,
\label{TBT}
\end{align}
where $n$ is the index of the unit cell,  and   $a^{}_{n}=\{a^{}_{n\Uparrow},a^{}_{n\Downarrow}\}$ and $b^{}_{n}=\{b^{}_{n\Uparrow},b^{}_{n\Downarrow}\}$  are spinors.

For a spinless DQD chain, the tunneling matrices in Eqs.~(\ref{T-IN}) and (\ref{T-EX}) become real numbers, $t^{}_{{\rm in},0}$ and $t^{}_{{\rm ex},0}$, and the effective tight-binding Hamiltonian  is identical to that of the single-level Su-Schrieffer-Heeger model~\cite{Su1979}
\begin{align}
H^{}_{{\rm T},0}=\sum^{}_{n}\left(t^{}_{{\rm in},0 } a^{\dagger}_{n}b^{}_{n}+t^{}_{{\rm ex},0 } a^{\dagger}_{n+1}b^{}_{n}+{\rm h.c.}\right),
\label{TB0}
\end{align}
with $a^{}_{n}$ and $b^{}_{n}$ corresponding to the operators  for the lowest localized orbitals on the QD1 and QD2, respectively.

\section{The Bloch spectrum of the double-quantum-dot chain ~\label{S-2}}

\begin{figure}
\centering
\includegraphics[width=0.4\textwidth]{Figs2.pdf}
\caption{  The spatial distribution of the confinement potential $U(x)$ in a unit cell; separate the unit cell into four regions (denoted by $I$, $II$, $III$ and $IV$)  in which the Bloch functions are derived (see text).}
\label{Figs2}
\end{figure}

Here the detailed calculation of the Bloch bands of our system is given.
For a one-dimensional lattice obeying periodic boundary conditions, the eigenstates are  the Bloch wave functions
\begin{align}
\Psi^{}_{k}(x)=\exp\left(ikx\right)u^{}_{k}(x),
\end{align}
 where $k$ is the wave vector and $u^{}_{k}(x)$ is the cell-periodic Bloch wave function. Obviously, the corresponding energy spectrum $\varepsilon(k)$ is given by
 \begin{align}
 H \Psi^{}_{k}(x)=\varepsilon(k)\Psi^{}_{k}(x),
 \label{B-F}
 \end{align}
where
 \begin{align}
  H =\frac{p^{2}}{2m_{e}}-V^{}_{\rm c}+U(x)+\alpha p \sigma_{y}+\frac{\Delta^{}_{\rm z}}{2}\sigma^{}_{n}.
\end{align}

 Similarly to the derivation in Sec.~\ref{S-1}, the solutions of Eq.~(\ref{B-F}) are obtained from the eigenfunctions of the Schr\"{o}dinger equations with different gate potentials
\begin{align}
\widetilde{H}^{}_{\zeta=0,1,2}\Psi^{}_{k,\zeta }(x)=E \Psi_{k,\zeta}(x).
\label{BSZ}
\end{align}
where $\Psi^{}_{k,\zeta}(x)=\exp\left(ikx\right)u^{}_{k,\zeta}(x)$  and $\varepsilon(k)=E-V^{}_{\rm c}$ (measured from the chemical potential $V_{\rm c}$).
 Then, analogously to the eigenstates in Eq.~(\ref{LQD}),  there exist four-fold degenerate solutions to Eq.~(\ref{BSZ}) for a fixed $k$, of the form
 \begin{align}
u^{}_{k,\zeta,\lambda}(x)=\exp\big[i(\gamma^{}_{\zeta,\lambda}-k)x\big]\begin{pmatrix}
w^{}_{\zeta,\lambda}\\
d^{}_{\zeta,\lambda}
\end{pmatrix}
\label{BF}
\end{align}
with the wave vector $\gamma^{}_{\zeta,\lambda}$ and the two spin components ($u^{}_{\zeta,\lambda}$ and $d^{}_{\zeta,\lambda}$) determined by Eqs.~(\ref{wk}) and (\ref{wkc}). It is convenient to divide the unit cell into four parts,  as shown in Fig.~\ref{Figs2}(a). The Bloch functions in each segment are expanded as
\begin{align}
u^{}_{k,I} (x)=&\sum^{4}_{\lambda=1}c^{}_{I,\lambda}u^{}_{k,0,\lambda} (x)~~~~ u^{}_{k,II} (x)=\sum^{4}_{\lambda=1}c^{}_{II,\lambda}u^{}_{k,1,\lambda}(x) \nonumber\\
u^{}_{k,III} (x)=&\sum^{4}_{\lambda=1}c^{}_{III,\lambda}u^{}_{k,0,\lambda}(x) ~~~~ u^{}_{k,IV} (x)=\sum^{4}_{\lambda=1}c^{}_{IV,\lambda}u^{}_{k,2,\lambda}(x),
\label{4u}
\end{align}
with $c^{}_{\Lambda,\lambda}$ $(\Lambda=I,~II,~III,~IV)$ representing the, yet to be determined, sixteen coefficients. The entire Bloch wave function is then
\begin{align}
u^{}_{k}(x)=&\Theta(x+L-d)[1-\Theta(x+d)]u^{}_{k,I}(x)+[1-\Theta(|x|-d)]u^{}_{k,II}(x)+[1-\Theta(x-L+d)]\nonumber\\
&\times  \Theta(x-d) u^{}_{k,III}(x) + \Theta(|x|-L+d)[1-\Theta(|x|-L)]u^{}_{k,IV}(x).
\label{bfc}
\end{align}
Combined with the periodicity of the cell Bloch functions, i.e., $u^{}_{k}(x_{2})=u^{}_{k}(x_{5})$ and $u^{\prime}_{k}(x_{2})=u^{\prime}_{k}(x_{5})$ with $x^{}_{5}\equiv x^{}_{2}+2L= L+d$, the energy bands are determined by the continuity of the wave function at the interfaces, $\lim_{\epsilon \rightarrow 0 }[u^{}_{k}(\xi+\epsilon)-u^{}_{k}(\xi+\epsilon)]=0$ and $\lim_{\epsilon \rightarrow 0 }[u^{\prime}_{k}(\xi+\epsilon)-u^{\prime}_{k}(\xi+\epsilon)]=0$,  with $\xi=x^{}_{1},x^{}_{3},x^{}_{4}$. Regarding the coefficients of the combinations in Eqs. (\ref{4u}) as a variable vector, the prerequisite conditions can be written in a matrix form
\begin{align}
\mathbf{W}(k,E) \cdot \hat{s}=0,
\label{MAT2}
\end{align}
where the vector is
$\hat{s}^{T}=\{c^{}_{I,1},c^{}_{I,2},c^{}_{I,3},c^{}_{I,4},c^{}_{II,1},c^{}_{II,2},c^{}_{II,3},c^{}_{II,4},c^{}_{III,1},c^{}_{III,2},c^{}_{III,3},c^{}_{III,4},c^{}_{IV,1},c^{}_{IV,2},c^{}_{IV,3},c^{}_{IV,4}\}$ and the matrix is
\begin{footnotesize}
\begin{align}
\mathbf{W}=\left(\begin{array}{cccccccccccccccc}
u^{(1,\uparrow)}_{k,0,\lambda}&u^{(1,\uparrow)}_{k,0,\lambda}&u^{(1,\uparrow)}_{k,0,\lambda}&u^{(1,\uparrow)}_{k,0,\lambda}&-u^{(1,\uparrow)}_{k,1,\lambda}&-u^{(1,\uparrow)}_{k,1,\lambda}&-u^{(1,\uparrow)}_{k,1,\lambda} &-u^{(1,\uparrow)}_{k,1,\lambda} &0&0&0&0 &0&0 &0&0\\
u^{(1,\downarrow)}_{k,0,\lambda}&u^{(1,\downarrow)}_{k,0,\lambda}&u^{(1,\downarrow)}_{k,0,\lambda}&u^{(1,\downarrow)}_{k,0,\lambda}&-u^{(1,\downarrow)}_{k,1,\lambda}&-u^{(1,\downarrow)}_{k,1,\lambda}&-u^{(1,\downarrow)}_{k,1,\lambda} &-u^{(1,\downarrow)}_{k,1,\lambda} &0&0&0&0 &0&0 &0&0\\
0&0&0&0&-u^{(3,\uparrow)}_{k,1,\lambda}&-u^{(3,\uparrow)}_{k,1,\lambda}&-u^{(3,\uparrow)}_{k,1,\lambda} &-u^{(3,\uparrow)}_{k,1,\lambda} &u^{(3,\uparrow)}_{k,0,\lambda}&u^{(3,\uparrow)}_{k,0,\lambda}&u^{(3,\uparrow)}_{k,0,\lambda}&u^{(3,\uparrow)}_{k,0,\lambda} &0&0 &0&0\\
0&0&0&0&-u^{(3,\downarrow)}_{k,1,\lambda}&-u^{(3,\downarrow)}_{k,1,\lambda}&-u^{(3,\downarrow)}_{k,1,\lambda} &-u^{(3,\downarrow)}_{k,1,\lambda} &u^{(3,\downarrow)}_{k,0,\lambda}&u^{(3,\downarrow)}_{k,0,\lambda}&u^{(3,\downarrow)}_{k,0,\lambda}&u^{(3,\downarrow)}_{k,0,\lambda} &0&0 &0&0\\
0&0&0&0&0&0&0 &0&u^{(4,\uparrow)}_{k,0,\lambda}&u^{(4,\uparrow)}_{k,0,\lambda}&u^{(4,\uparrow)}_{k,0,\lambda}&u^{(4,\uparrow)}_{k,0,\lambda} &-u^{(4,\uparrow)}_{k,2,\lambda}&-u^{(4,\uparrow)}_{k,2,\lambda}&-u^{(4,\uparrow)}_{k,2,\lambda} &-u^{(4,\uparrow)}_{k,2,\lambda}\\
0&0&0&0&0&0&0 &0&u^{(4,\downarrow)}_{k,0,\lambda}&u^{(4,\downarrow)}_{k,0,\lambda}&u^{(4,\downarrow)}_{k,0,\lambda}&u^{(4,\downarrow)}_{k,0,\lambda} &-u^{(4,\downarrow)}_{k,2,\lambda}&-u^{(4,\downarrow)}_{k,2,\lambda}&-u^{(4,\downarrow)}_{k,2,\lambda} &-u^{(4,\downarrow)}_{k,2,\lambda}\\
u^{(2,\uparrow)}_{k,0,\lambda}&u^{(2,\uparrow)}_{k,0,\lambda}&u^{(2,\uparrow)}_{k,0,\lambda}&u^{(2,\uparrow)}_{k,0,\lambda}&0&0&0 &0 &0&0&0&0 &-u^{(5,\uparrow)}_{k,2,\lambda}&-u^{(5,\uparrow)}_{k,2,\lambda}&-u^{(5,\uparrow)}_{k,2,\lambda} &-u^{(5,\uparrow)}_{k,2,\lambda}\\
u^{(2,\downarrow)}_{k,0,\lambda}&u^{(2,\downarrow)}_{k,0,\lambda}&u^{(2,\downarrow)}_{k,0,\lambda}&u^{(2,\downarrow)}_{k,0,\lambda}&0&0&0 &0 &0&0&0&0 &-u^{(5,\downarrow)}_{k,2,\lambda}&-u^{(5,\downarrow)}_{k,2,\lambda}&-u^{(5,\downarrow)}_{k,2,\lambda} &-u^{(5,\downarrow)}_{k,2,\lambda}\\
\tilde{u}^{(1,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(1,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(1,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(1,\uparrow)}_{k,0,\lambda}&-\tilde{u}^{(1,\uparrow)}_{k,1,\lambda}&-\tilde{u}^{(1,\uparrow)}_{k,1,\lambda}&-\tilde{u}^{(1,\uparrow)}_{k,1,\lambda} &-\tilde{u}^{(1,\uparrow)}_{k,1,\lambda} &0&0&0&0 &0&0 &0&0\\
\tilde{u}^{(1,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(1,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(1,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(1,\downarrow)}_{k,0,\lambda}&-\tilde{u}^{(1,\downarrow)}_{k,1,\lambda}&-\tilde{u}^{(1,\downarrow)}_{k,1,\lambda}&-\tilde{u}^{(1,\downarrow)}_{k,1,\lambda} &-\tilde{u}^{(1,\downarrow)}_{k,1,\lambda} &0&0&0&0 &0&0 &0&0\\
0&0&0&0&-\tilde{u}^{(3,\uparrow)}_{k,1,\lambda}&-\tilde{u}^{(3,\uparrow)}_{k,1,\lambda}&-\tilde{u}^{(3,\uparrow)}_{k,1,\lambda} &-\tilde{u}^{(3,\uparrow)}_{k,1,\lambda} &\tilde{u}^{(3,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(3,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(3,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(3,\uparrow)}_{k,0,\lambda} &0&0 &0&0\\
0&0&0&0&-\tilde{u}^{(3,\downarrow)}_{k,1,\lambda}&-\tilde{u}^{(3,\downarrow)}_{k,1,\lambda}&-\tilde{u}^{(3,\downarrow)}_{k,1,\lambda} &-\tilde{u}^{(3,\downarrow)}_{k,1,\lambda} &\tilde{u}^{(3,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(3,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(3,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(3,\downarrow)}_{k,0,\lambda} &0&0 &0&0\\
0&0&0&0&0&0&0 &0&\tilde{u}^{(4,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(4,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(4,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(4,\uparrow)}_{k,0,\lambda} &-\tilde{u}^{(4,\uparrow)}_{k,2,\lambda}&-\tilde{u}^{(4,\uparrow)}_{k,2,\lambda}&-\tilde{u}^{(4,\uparrow)}_{k,2,\lambda} &-\tilde{u}^{(4,\uparrow)}_{k,2,\lambda}\\
0&0&0&0&0&0&0 &0&\tilde{u}^{(4,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(4,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(4,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(4,\downarrow)}_{k,0,\lambda} &-\tilde{u}^{(4,\downarrow)}_{k,2,\lambda}&-\tilde{u}^{(4,\downarrow)}_{k,2,\lambda}&-\tilde{u}^{(4,\downarrow)}_{k,2,\lambda} &-\tilde{u}^{(4,\downarrow)}_{k,2,\lambda}\\
\tilde{u}^{(2,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(2,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(2,\uparrow)}_{k,0,\lambda}&\tilde{u}^{(2,\uparrow)}_{k,0,\lambda}&0&0&0 &0 &0&0&0&0 &-\tilde{u}^{(5,\uparrow)}_{k,2,\lambda}&-\tilde{u}^{(5,\uparrow)}_{k,2,\lambda}&-\tilde{u}^{(5,\uparrow)}_{k,2,\lambda} &-\tilde{u}^{(5,\uparrow)}_{k,2,\lambda}\\
\tilde{u}^{(2,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(2,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(2,\downarrow)}_{k,0,\lambda}&\tilde{u}^{(2,\downarrow)}_{k,0,\lambda}&0&0&0 &0 &0&0&0&0 &-\tilde{u}^{(5,\downarrow)}_{k,2,\lambda}&-\tilde{u}^{(5,\downarrow)}_{k,2,\lambda}&-\tilde{u}^{(5,\downarrow)}_{k,2,\lambda} &-\tilde{u}^{(5,\downarrow)}_{k,2,\lambda}
\end{array}\right)
\label{CMAT}
\end{align}
\end{footnotesize}
with $u^{(j,\uparrow)}_{k,\zeta,\lambda}= \exp[i(\gamma^{}_{\zeta,\lambda}-k)x_{j}]w^{}_{\zeta,\lambda}$ ($j=1,2,3,4,5$), $u^{(j,\downarrow)}_{k,\zeta,\lambda}= \exp[i(\gamma^{}_{\zeta,\lambda}-k)x_{j}]d^{}_{\zeta,\lambda}$, $\tilde{u}^{(j,\uparrow)}_{k,\zeta,\lambda}= (\gamma^{}_{\zeta,\lambda}-k)\exp[i(\gamma^{}_{\zeta,\lambda}-k)x_{j}]w^{}_{\zeta,\lambda}$, and $\tilde{u}^{(j,\downarrow)}_{k,\zeta,\lambda}= (\gamma^{}_{\zeta,\lambda}-k)\exp[i(\gamma^{}_{\zeta,\lambda}-k)x_{j}]d^{}_{\zeta,\lambda}$. A nontrivial solution to Eq.~(\ref{MAT2}) requires
\begin{align}
\det [\mathbf{W} (k,E)]=0,
\label{GF}
\end{align}
 and yields an implicit equation for $k$ and $E$ of Eqs.~(\ref{wk}),~(\ref{wkc}), and ~(\ref{CMAT}). The exact energy spectrum of the Bloch band $\varepsilon (k)=E-V^{}_{\rm c}$ is found by solving numerically for $E$ at a fixed value of $k$. Substituting the so-derived values of $E$ in Eq.~(\ref{MAT2}) yields the values of the coefficients $c^{}_{\Lambda,\lambda}$, i.e., gives the specific form of the Bloch wave functions, with the normalization condition $\int^{L}_{-L}u^{\ast}_{k}(x)u^{ }_{k}(x)dx=1$.

In the absence of a magnetic field or when $\theta=0$, the two spin degrees of freedom are separable and the number of the corresponding coefficients  for ascertaining the Bloch band is reduced by half, which is similar to a spinless chain.  In this case, the matrix $\mathbf{W}(k,E)$ is  decomposed into two independent submatrices whose determinants yield the Bloch spectrum, as shown in Fig.3 (a) for the case of no magnetic field.

\section{Edge states of a finite double-quantum-dot chain}
\label{S-3}

\begin{figure}
\centering
\includegraphics[width=0.68\textwidth]{Figs3.pdf}
\caption{  The energy  spectrum (as a function of $k$) of a double quantum-dot chain made of  InAs, for different values of the magnetic field and in the absence of the spin-orbit coupling. Here $\Delta V=10$~meV and $\theta=\pi/2$. (a)  Small magnetic fields,    $\Delta^{}_{\rm z}<t^{}_{\rm ex,0 }-t^{}_{\rm in,0 }$. (b) Large  magnetic fields,  $\Delta^{}_{\rm z}>t^{}_{\rm ex,0 }+t^{}_{\rm in,0 }$. $k$ is measured in units of $1/(2L)$ and the energies of two pure spin edge states, $\varepsilon^{}_{\uparrow}$ and $\varepsilon^{}_{\downarrow}$, are also indicated.}
\label{Figs3}
\end{figure}

Here we consider the edge modes of an open  chain of DQDs.
 We begin with the configuration where there is no Zeeman field,
and the electron spin is polarized solely by the field  induced by the spin-orbit interaction.
In this case the tight-binding Hamiltonian reads
 \begin{align}
H^{}_{\rm T,S}=\sum^{}_{n} (t^{}_{\rm in,0}a^{\dagger}_{n}e^{i\varphi^{}_{\rm in}\sigma^{}_{z}}b^{}_{n}+t^{}_{\rm ex,0}a^{\dagger}_{n+1}e^{i\varphi^{}_{\rm ex}\sigma^{}_{z}}b^{}_{n}+{\rm h.c.} ),
\nonumber
 \end{align}
and there are no spin-flip transitions.

The energy eigenfunctions in an  open DQD chain are distinguished by the different quasi-spin states
\begin{align}
|\Psi^{}_{\tau=\Uparrow,~\Downarrow }\rangle=\sum^{N}_{n=1}\left( A^{}_{n,\tau} a^{\dagger}_{n\tau}+B^{}_{n,\tau}b^{\dagger}_{n\tau}\right)|0\rangle,
\label{F-LS}
\end{align}
where $|0\rangle$ is the vacuum state.  The coefficients
$A^{}_{n,\tau}$ and $B^{}_{n,\tau}$
 are determined by the spin-dependent Schr\"{o}dinger equation, $H^{}_{\rm T,S}|\Psi^{}_{\tau}\rangle=E|\Psi^{}_{\tau}\rangle$, where $E^{}_{}$ is  the energy eigenvalue. Specifically, the Schr\"{o}dinger equation for the quasi-spin up state can be written as
\begin{align}
\sum^{}_{n}\left( t^{}_{\rm in, 0}e^{i\varphi^{}_{\rm in}}a^{\dagger}_{n\Uparrow} b^{}_{n\Uparrow}+ t^{}_{\rm ex,0}e^{i\varphi^{}_{\rm ex}}a^{\dagger}_{n+1\Uparrow}b^{}_{n\Uparrow}+{\rm h.c.}\right)|\Psi^{}_{ \Uparrow }\rangle=E^{}_{ }|\Psi^{}_{ \Uparrow}\rangle.
\label{EGS}
\end{align}
Then, the explicit Shr\"{o}dinger equations in  the $n$-th cell are
\begin{align}
 t^{}_{\rm in, 0}e^{i\varphi^{}_{\rm in}}B^{}_{n,\Uparrow}+  t^{}_{\rm ex, 0}e^{i\varphi^{}_{\rm ex}}B^{}_{n-1,\Uparrow}=E^{}_{} A^{}_{n,\Uparrow}\ ,~~~~~~~~~ t^{}_{\rm in, 0}e^{-i\varphi^{}_{\rm in}}A^{}_{n,\Uparrow}+ t^{}_{\rm ex, 0}e^{-i\varphi^{}_{\rm in}} A^{}_{n+1,\Uparrow} =E^{}_{}B^{}_{n,\Uparrow}\ .
 \label{CQES}
 \end{align}
These can be reduced to a recursion formula when $E^{}_{ } \neq 0$
\begin{align}
X^{}_{n,\Uparrow}\left(t^{2}_{\rm in,0}+t^{2}_{\rm ex, 0}-E_{}^{2}\right)+ t^{}_{\rm in,0}t^{}_{\rm ex,0}e^{i(\varphi^{}_{\rm in}-\varphi^{}_{\rm ex})}X^{}_{n+1,\Uparrow}+t^{}_{\rm in,0}t^{}_{\rm ex,0}e^{i(\varphi^{}_{\rm ex}-\varphi^{}_{\rm in})}X^{}_{n-1,\Uparrow}=0\ ,
\end{align}
with $X=A,~B$. This recursion relation is solved by  the transfer-matrix scheme
\begin{align}
\begin{pmatrix}
X_{n+1,\Uparrow}\\
X_{n,\Uparrow}
\end{pmatrix}=\mathbf{T}\begin{pmatrix}
X_{n,\Uparrow}\\
X_{n-1,\Uparrow}
\end{pmatrix}\ ,
\end{align}
with
\begin{align}
\mathbf{T}=\begin{pmatrix}
 m^{}_{11}& m^{}_{12}\\
1&0
\end{pmatrix}\ ,
\end{align}
in which $m^{}_{11} =e^{i(\varphi^{}_{\rm ex}-\varphi^{}_{\rm in})}(E^{2}_{}-t^{2}_{\rm in,0}-t^{2}_{\rm ex, 0})/(t^{}_{\rm in,0}t^{}_{\rm ex,0})$ and $m^{}_{12}=e^{2i(\varphi^{}_{\rm ex}-\varphi^{}_{\rm in})}$. The formation of edge modes hinges on the eigenvalues $\lambda^{}_{\pm}= \frac{1}{2} \left(m^{}_{11}\pm\sqrt{m^{2}_{11}+4m^{}_{12}}\right)$ of the transfer matrix $\mathbf{T}$. An eigenvalue with an absolute value larger than 1, $|\lambda^{}_{+}|>1$, excludes the possibility to form edge modes with nonzero energy~\cite{Li2015}.
Thus, the edge modes only occur for $E^{}_{ }=0$. The recursion relations for this case are
\begin{align}
t^{}_{\rm in, 0}e^{-i\varphi^{}_{\rm in }} A_{n,\Uparrow}+t^{}_{\rm ex,0} e^{-i\varphi^{}_{\rm ex }} A_{n+1,\Uparrow}=0\ ,~~~~~~~~  t_{\rm in, 0} e^{i\varphi^{}_{\rm in }}B_{n,\Uparrow}+t_{\rm ex, 0}e^{i\varphi^{}_{\rm ex }} B_{n-1,\Uparrow}=0 .
\end{align}
As occurs for the Su-Schrieffer-Heeger model~\cite{Su1979}, the zero-energy edge modes appear only when $ t^{}_{\rm in,0}/t^{}_{\rm ex, 0}<1$. By introducing the notation  $ r^{}_{ }=(t^{}_{\rm in,0}/t^{}_{\rm ex, 0})e^{i(\varphi^{}_{\rm in }-\varphi^{}_{\rm ex })}$, the zero-energy edge modes are given by
\begin{align}
|E^{}_{a,\Uparrow}\rangle=A^{}_{1,\Uparrow}\sum^{N}_{n=1}(r_{}^{\ast})^{n-1}a^{\dagger}_{n\Uparrow}|0\rangle\ ,  ~~~~~~~~ |E^{}_{b,\Uparrow}\rangle=B^{}_{N,\Uparrow}\sum^{N}_{n=1}r_{}^{ N-n}b^{\dagger}_{n\Uparrow}|0\rangle,
\label{SOE}
\end{align}
with the coefficients $A^{}_{1,\Uparrow}$ and $B^{}_{N,\Uparrow}$  determined by the normalization conditions  $\langle E^{}_{a/b,\Uparrow} | E^{}_{a/b ,\Uparrow}\rangle=1$. Meanwhile, in the case of a small $N$, the overlap between $|E^{}_{a,\Uparrow}\rangle$ and $|E^{}_{b,\Uparrow}\rangle$    results in the hybridization between the two zero-energy edge modes, as indicated in Fig.~3(b) for an open spinless chain with $N=20$.
Analogously, the quasi-spin down edge states can be derived by the same method as above, and all the zero-energy edge modes are protected by  the chiral symmetry as explained in the main text.

 The situation is changed when a Zeeman interaction is applied. We first consider the Hamiltonian in the absence of spin-orbit coupling, i.e., for
 $\varphi^{}_{\rm in/ex}=0$,
\begin{align}
 H^{ }_{  \rm  T,Z}= \sum^{}_{n}\frac{\Delta^{}_{\rm z}}{2}\left(a^{\dagger}_{n}\sigma^{}_{z}a^{}_{n}+b^{\dagger}_{n}\sigma^{}_{z}b^{}_{n}\right)+\sum^{}_{n}\sum^{}_{\tau=\uparrow,\downarrow}\left(t^{}_{\rm in,0} a^{\dagger}_{n\tau}b^{}_{n\tau}+t^{}_{\rm ex, 0}a^{\dagger}_{n+1\tau}b^{}_{n\tau}+ \rm{h.c.}\right).
 \label{TBP}
  \end{align}
  Notably, because there is no spin-orbit interaction, in this case $\tau=\uparrow,\downarrow$ represent the two pure spin states  polarized by the Zeeman field.
 In Fourier space, this   Hamiltonian reads
\begin{align}
H^{}_{\rm T,Z}(k)= \begin{pmatrix}
\Delta^{}_{\rm z}/2 &t^{}_{\rm  in, 0 }+e^{-ik}t^{}_{\rm  ex,0 }&0&0\\
t^{ }_{\rm  in, 0}+e^{ik}t^{ }_{\rm ex,0} &\Delta^{}_{\rm z}/2&0&0\\
0&0&-\Delta^{}_{\rm z}/2& t^{}_{\rm  in,0 }+e^{-ik}t^{}_{\rm  ex,0}\\
0&0&t^{}_{\rm  in,0 }+e^{ik}t^{\ast }_{\rm  ex,0} &-\Delta^{}_{\rm z}/2
\end{pmatrix},
\end{align}
and the energies of the spin-dependent Bloch bands are
\begin{align}
E^{ }_{\uparrow \pm} (k)=& \left(\Delta^{}_{\rm z}/2\right)\pm \sqrt{t^{2}_{\rm in,0}+ t^{2}_{\rm ex,0} +2t^{ }_{\rm in,0}t^{ }_{\rm ex,0}\cos k}\nonumber\\
E^{ }_{\downarrow \pm} (k)=& -\left(\Delta^{}_{\rm z}/2\right) \pm \sqrt{t^{2}_{\rm in,0}+ t^{2}_{\rm ex,0} +2t^{ }_{\rm in,0}t^{ }_{\rm ex,0}\cos k}.
\label{enes}
 \end{align}
 As seen, the Hamiltonian Eq.~(\ref{TBP}) separates  into two independent   models, each being described by   a Su-Schrieffer-Heeger model with a finite on-site energy. Consequently, one expects two pairs of  edge modes in the topological regime in which $t^{}_{\rm ex,0}>t^{}_{\rm in,0}$. However, when
\begin{align}
t^{}_{\rm ex,0}- t^{}_{\rm in,0}\leq\Delta^{} _{\rm z} \leq t^{}_{\rm ex,0}+ t^{}_{\rm in,0},
\label{BBC}
\end{align}
the overlap of an edge mode belonging to
a certain spin projection with the bulk of states of the opposite projection
masks its topological character.
It follows that in the absence of spin-orbit coupling the edge states reside within the energy gap provided that
 $ \Delta^{}_{\rm z}<t^{}_{\rm ex,0}-t^{}_{\rm in,0} $ at small magnetic fields, and  $ \Delta^{}_{\rm z}>t^{}_{\rm ex,0}+t^{}_{\rm in,0} $ at large  fields, as shown in Fig.~ \ref{Figs3}.

When both the spin-orbit coupling and the Zeeman magnetic field are active, the energy eigenfunctions of an open dimerized chain are  the quasi-spin mixed states
\begin{align}
|\Psi^{}_{}\rangle =\sum^{N}_{n=1} \left(A^{}_{n,\Uparrow}a^{\dagger}_{n\Uparrow}+ A^{}_{n,\Downarrow}a^{\dagger}_{n\Downarrow}+B^{}_{n,\Uparrow}b^{\dagger}_{n\Uparrow}+B^{}_{n,\Downarrow}b^{\dagger}_{n\Downarrow}\right)|0\rangle.
\label{PSF0}
\end{align}
The coefficients are determined
 by the Schr\"{o}dinger equation
$H^{}_{\rm T}|\Psi^{}_{}\rangle =E_{ }|\Psi^{}_{}\rangle$, with the discrete Hamiltonian $H^{}_{\rm T}$ given in Eq.~(\ref{TBT}).
The edge modes are distinguished by the spatial density distributions of the state, i.e., they are residing on either side of the chain. We quantify this property by choosing
 \begin{align}
 \eta^{}_{\rm L} =\sum^{3}_{n=1}P^{}_{n}+\sum^{N}_{n=N-2}P^{}_{n}
\end{align}
with $P^{}_{n}= |A^{}_{n,\Uparrow}|^{2}+|B^{}_{n,\Uparrow}|^{2}+ |A^{}_{n,\Downarrow}|^{2}+|B^{}_{n,\Downarrow}|^{2}$, $n=1,2,3$ indicating the leftmost three sites of the chain,  and $n=N-2,N-1,N$ corresponding to rightmost three sides.  For an open InAs chain of $N=20$ sites and $\Delta V=10$~meV, the specifical values of $\eta^{}_{\rm L}$ for the edge states under different magnetic fields are scaled by the color bar shown in Fig.~\ref{Figs4}. In order to count the number of  edge modes as more as possible, here  we assume the presence of an edge state requires $\eta^{}_{\rm L}\geq 0.45$, with the  threshold value determined by the probability density distribution of a predestined edge state in  the prototype SSH model with a small degree of dimerization, $\Delta V=1.0$~meV.

Being consistent with the results shown in Fig.~3(d), the numerical values of  $\eta^{}_{\rm L}$  in Fig~\ref{Figs4}(a) confirm  the appearance of  edge states  near the zero degrees of dimerization ( i.e., $\Delta V\rightarrow 0^{+}$), but characterized with the weak localization.   Under the interplay between the spin-orbit and Zeeman interactions, the large locality  parameter $\eta^{}_{\rm L}$  in Fig.~\ref{Figs4}(b) demonstrates the presence of additional edge modes living in the gap induced by the anti-crossing of the energy levels. However, in contrast to the Zeeman-split edge states, these  extra edge modes are not stable and can be merged into the bulk states by decreasing the tilting angle $\theta$ for $0\leq\theta<\pi/2$, and that is  manifested by the reduction of $\eta^{}_{\rm L}$ shown in Fig.~\ref{Figs4}(c). Besides,  the insensitivity of the energy values to $\Delta^{}_{\rm z}$, seen in Fig~\ref{Figs4}(b), can also be illustrated by  the high degree of spin mixing intrinsic to the additional edge states, as reflected by the small spin  density distribution, i.e., $S_{n}=\left(|A^{}_{n\Uparrow}|^{2}+|B^{}_{n\Uparrow}|^{2}-|A^{}_{n\Downarrow}|^{2}-|B^{}_{n\Downarrow}|^{2}\right)/2$, shown in Fig.~\ref{Figs4}(d).

\begin{figure}
\centering
\includegraphics[width=0.7\textwidth]{Figs4.pdf}
\caption{ (a) The energy spectrum of an open InAs  chain comprising 20 sites as a function of the potential difference $\Delta V$, with the magnetic field $B=0.2$~T, the tilting angle $\theta=0.5\pi$, and the spin-orbit length $x^{}_{\rm so}=180$~nm. Moreover, for each edge state, the locality parameter $\eta^{}_{\rm L}$ is scaled by the color bar shown in the right side. (b) The energy spectrum of the open InAs chain as a function of the Zeeman splitting $\Delta_{\rm z}$, with $\Delta V=10$~nm, $\theta=0.5\pi$, and  $x^{}_{\rm so}=180$~nm. (c) For the case of $B=1.0$~T, the energy spectrum of the open chain as a function of $\theta$, when $\Delta V=10$~nm and $x^{}_{\rm so}=180$~nm.  Moreover, the spatial spin density distributions  of the two  edge modes, as indicated by a square and a diamond in panel (c), are displayed in (d).  }
\label{Figs4}
\end{figure}

 \footnotetext[1]{asb}