\section{Introduction}

Object matching techniques have been widely used in many fields of computer vision, including 2D and 3D image analysis, object recognition, biomedical identification, and object tracking. There is a rich literature on finding meaningful approximate isomorphisms between pair of objects that are represented as graphs~\cite{umeyama88,cour2006,emms2009,gori2005,wilson2008,egozi2012,Hu2013,Hu2014}. Many tasks, however, require to solve the so-called multi-object matching problem, i.e., finding consistent maps among all pairs of objects within a collection. Examples include non-rigid structure from motion \cite{Bregler2000, Dai2012} and shared object discovery \cite{ChenXL2014}. In this context, a central task is how to utilize the data collection as a regularizer to improve the maps computed between pairs of objects in isolation~\cite{Huang2013,Chen2014, Zhou2015}.

A generic constraint that one can utilize to improve maps among a collection is the so-called \emph{cycle consistency} constraint, namely composition of maps along any two paths sharing the same starting and end objects are identical. A technical challenge of utilizing this constraint is that it is impossible to check all cycles for consistency, due to the fact that the number of paths increase exponentially with the total number of objects. Recent works on joint matching have shown that the cycle consistency constraint can be translated into a much more manageable constraint, i.e., the data matrix that stores pair-wise maps in blocks is positive semidefinite and low-rank \cite{Huang2013,Pachauri2013}. Based on this connection, people have formulated multi-object matching as solving semidefinite programs (or SDP), which are convex relaxations of the corresponding matrix recovery problem. These algorithms achieved near-optimal exact recovery conditions~\cite{Huang2013, Chen2014}. On the other hand, solving semidefinite programs are computationally expensive. In a recent work, Zhou et al~\cite{Zhou2015} attempt to address the computational issue using alternating minimization for efficient low-rank matrix recovery.

In this paper, we propose a novel framework that utilizes the cycle-consistency constraint in a hierarchical manner for scalable multiple object matching. We show how to apply this framework to extend the methods described in~\cite{Huang2013,Chen2014,Zhou2015}.
In particular, instead of jointly imposing the global consistency constraint among all pair-wise maps~\cite{Huang2013, Chen2014, Zhou2015}, we split the input object collection into overlapping subsets, and impose consistency within each subset. We then impose consistency between maps across the subsets. Interestingly, we show that by combing these two consistency constraints together, we can guarantee global consistency under mild conditions (See Section \ref{sec:consistency}). Yet computationally, such a decoupled approach yields significant performance gains, when compared with existing approaches.

\subsection{Related Work}

Early works on multi-object matching (e.g.,~\cite{zass2008, Lee2011}) extend pairwise matching schemes to the multi-object setting without explicitly considering the map consistency constraint. \cite{Zach2010, Nguyen2011} proposed to detect inconsistent cycles, and formulate multi-object matching as solving combinatorial optimizations, i.e., removing bad maps to break all inconsistent cycles. Recently, people have proposed to formulate non-convex optimization problems by using the cycle
consistency constraint as an explicit constraint for either pixel-wise flow computation \cite{Zhou2015}, sparse feature matching \cite{Yan2015}, sparse shape modeling \cite{Cosmo2016}, or structure from motion \cite{Pachauri2014}. These problems are, as a consequence, hard to solve and do not admit exact recovery conditions. Recent works
\cite{Huang2013, Pachauri2013} showed that consistent maps could be extracted
from the spectrum of a data matrix that encodes pair-wise maps in blocks. Along this line of research, Huang and Guibas \cite{Huang2013} proposed an elegant solution by formulating the problem as convex relaxation and discussed the theoretical conditions for exact recovery. The result is further analyzed in \cite{Chen2014} under the condition that the underlying rank of the variable matrix is known or can be reliably estimated. Yan et al.~\cite{Yan_ICCV2015, Yan_TIP2015} also proposed matrix factorization based methods to enforce the cycle-consistency constraint. These methods, however, are not scalable to large-scale datasets, due to the cost of solving semidefinite programs. Zhou et al.~\cite{Zhou2015} enforce the positive semidefinite constraint using explicit low-rank factorizations, leading to improved computational efficiency. In contrast to these methods, our approach opens a new direction to  enforcing the cycle-consistency constraint, i.e., by splitting the datasets into overlapping subsets. This leads to further improvements in terms of computational efficiency. Most recently, Leonardos et. al. \cite{Leonardox2017} proposed a distributed consensus-based algorithm as an extension of \cite{Pachauri2013}. Their method, however, cannot handle partial matches.



We organize the reminder of this paper as follows. First, we discuss the problem setup and analyze the conditions in Section~\ref{sec:consistency}. Second, we discuss the formulation of our approach in Section~\ref{sec:formulation}. In Section~\ref{sec:solver}, we present an alternating direction method of multipliers (ADMM) for solving the induced optimization problem, leading to a parallel algorithm via generalized message passing. Last but not the least, we demonstrate the effectiveness of our approach on both synthetic and real examples in Section~\ref{sec:exp}.


\section{Consistency}
\label{sec:consistency}

In this section, we extend the cycle-consistency formulation described in~\cite{Huang2013} to the distributed setting. The key result is a sufficient condition on which cycle-consistency among sub-collections induces global cycle-consistency.

We begin with introducing the notations that are necessary to formally state this sufficient condition. For simplicity, we assume maps between objects are given by permutations. However, the argument can be easily extended to the case where objects are partially similar with each other. Formally speaking, we consider a map graph ${\cal G} = ({\cal V} = (H_1, \cdots, H_{n}), {\cal E})$. The vertex set ${\cal V}$ consists of objects to be matched, and each object $H_i$ is given by $m$ points (e.g., key points extracted from an image). The edge set ${\cal E}$ connects a subset of pairs of objects. Each edge $(i,j)\in{\cal E}$ is associated with a permutation $\phi_{ij}: H_i \rightarrow H_j$. We first define the global consistency of $\phi_{ij}, \forall (i,j)\in {\cal E}$:





\begin{mydef}[Cycle Consistency]
A map graph ${\cal G}=\{{\cal V},{\cal E}\}$ is \emph{cycle consistent} if for every node $v_i$ and every cycle $v_i-v_{i_1}-\cdots-v_{i_k}-v_i$, the composite map along this cycle is the identity map, i.e.
\[
\phi_{ii_1}\circ\cdots\circ\phi_{i_k i} = \emph{identity}.
\]
\end{mydef}

Now we introduce an equivalent formulation of enforcing the cycle-consistency among ${\cal G}$ as enforcing the cycle-consistency among subgraphs of ${\cal G}$, if these subgraphs satisfy certain conditions. Towards this end, we introduce two conditions among collection of subgraphs of ${\cal G}$. The first condition concerns a pair of sub-graphs:

\begin{mydef}[Joint Normal]
Let ${\cal G}_i=\{{\cal V}_i,{\cal E}_i\}$,  ${\cal G}_j=\{{\cal V}_j,{\cal E}_j\}$ be the two subgraphs of ${\cal G} = \{{\cal V}, {\cal E}\}$. We say ${\cal G}_i$ and ${\cal G}_j$ are \textsl{joint normal} if the vertex sets ${\cal V}_i \setminus {\cal V}_j$ and ${\cal V}_j \setminus {\cal V}_i$ are not connected by any edge of ${\cal E}$:
\[
(s,t)\notin{\cal E}, ~\forall H_s\in{\cal V}_i\backslash{\cal V}_j, H_t\in{\cal V}_j\backslash{\cal V}_i
\]
\end{mydef}

\begin{figure}
\begin{subfigure}{.48\linewidth} \centering \includegraphics[width=0.96\linewidth]{joint_normal}
\caption{{\it joint normal}}
\end{subfigure}
\begin{subfigure}{.48\linewidth} \centering \includegraphics[width=0.96\linewidth]{joint_normal_counter}
\caption{Counter {\it joint normal}}
\end{subfigure}
\caption{Examples of subgraphs that are (a) joint normal and (b) not joint normal.}
\label{fig:joint_normal}
\end{figure}
As illustrated in Figure~\ref{fig:joint_normal}, two subgraphs ${\cal G}_i=\{{\cal V}_i,{\cal E}_i\}$ and ${\cal G}_j=\{{\cal V}_j,{\cal E}_j\}$ are {\it joint normal} if their common subgraph $\{{\cal V}_i\cap {\cal V}_j,{\cal E}_i \cap {\cal E}_j\}$ is either (i) empty, or (ii) connected, and there is no edge between a vertex in one subgraph to a vertex in the other subgraph except those in the common subgraph. In contrast, the two subgraphs illustrated in Figure~\ref{fig:joint_normal}(b) are not joint normal since there exists an edge that connects the non-overlapping sets of these two subgraphs.

The second condition concerns a topological constraint among all the sub-graphs. We state this condition using the notation of simplicial complex as detailed below:
\begin{mydef}[Cover Complex]
Let $\{{\cal G}_i = ({\cal V}_i, {\cal E}_i), 1\leq i \leq K\}$ be a set of sub-graphs that cover ${\cal G}$, i.e., $\cup_{i=1}^{K}{\cal V}_i = {\cal V}$. We define the cover complex ${\cal K}$ of these sub-graphs $\{{\cal G}_i = ({\cal V}_i, {\cal E}_i), 1\leq i \leq K\}$ so that ${\cal K}$ collects every subset $ \{i_1, \cdots, i_{k}\}\subset \{1, \cdots, K\}$ if the intersections of ${\cal V}_{i_1}, \cdots, {\cal V}_{i_k}$ is non-empty, i.e. $\cap_j{\cal V}_{i_j}\neq\varnothing$.
\end{mydef}
 
We now state the decoupled cycle-consistency theorem that relates the global cycle consistency and the cycle consistency on each subgraph:
 
\begin{mythm}[Decoupled Cycle-Consistency]
\label{thm:consistency}
Let ${\cal G}$ be a map graph, ${\cal G}_1, \cdots, {\cal G}_{K}$ be a cover of ${\cal G}$, and ${\cal K}$ be the cover complex. Then ${\cal G}$ is cycle consistent if
\begin{enumerate}\denselist
\item ${\cal G}_i$ is \emph{cycle consistent} $\forall i$,
\item ${\cal G}_i$ and ${\cal G}_j$ are \emph{joint normal} $\forall (i,j)\in{\cal E}$,
\item ${\cal K}$ is \emph{simply connected} (c.f. \cite{Hatcher2002}).
\end{enumerate}
\end{mythm}

Here we say the cover complex ${\cal K}$ is \emph{simply connected} if every closed curve can be deformed to a point (or in other words the region enclosed by this curve has no-holes). Please refer to~\cite{Hatcher2002} for a more general definition. This theorem states that the cycle consistency property on each sub-graph would be propagated to the global consistency, if the cover complex ${\cal K}$ satisfies the conditions stated in Theorem~\ref{thm:consistency}.
The proof to Theorem \ref{thm:consistency} is left to the supplementary material.

Note that the $3^\textnormal{nd}$ condition in Theorem \ref{thm:consistency} is necessary. Figure~\ref{fig:cover_examples}(a) provides a simple counter example, which satisfies the $1^\textnormal{st}$ and $2^\textnormal{nd}$ conditions in Theorem \ref{thm:consistency}. The cover complex ${\cal K}$, however, is homologous to the Torus $T^2$, which is not simply connected. It is easy to see that the local-consistency (which is trivial as each sub-graph is given by an edge) does not lead to the global consistency among these three edges.

\begin{figure}
\begin{subfigure}{.96\linewidth}
\centering
\includegraphics[height=3cm]{graphs_triangle_cover}
\hspace{20pt}
\includegraphics[height=3cm]{triangle_empty}
\caption{Subgraphs that forms an empty triangle cover complex}
\end{subfigure}
\begin{subfigure}{.96\linewidth}
\centering
\includegraphics[height=3cm]{graphs_solid_triangle_cover}
\hspace{20pt}
\centering \includegraphics[height=3cm]{triangle_shape}
\caption{Subgraphs that form a solid triangle cover complex}
\end{subfigure}
\begin{subfigure}{.96\linewidth}
\centering
\includegraphics[height=3cm]{graphs_tetrahedron_cover}
\hspace{20pt}
\centering \includegraphics[height=3cm]{euclid_tetrahedron}
\caption{Subgraphs that form an empty tetrahedron cover complex}
\end{subfigure}
\caption{Given local consistency, (a) is not globally consistent, while (b) and (c) are guaranteed to be globally consistent. (Left: the map graph with subgraphs circled out. Right: the corresponding cover complex.)}
\label{fig:cover_examples}
\end{figure}

Figure~\ref{fig:cover_examples}(c) provides another example to understand the correctness of Theorem~\ref{thm:consistency}. In this case, there are four objects. It is clear that enforcing the cycle-consistency among all four triple sub-graphs induces the cycle-consistency on the original graph. This argument aligns with Theorem~\ref{thm:consistency} as ${\cal K}$ is simply connected. We defer detailed explanations to the supplementary material.

\begin{algorithm}[ht]
\caption{Greedy Construction of ${\cal K}$}
\label{alg:greedyK}
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{Map graph ${\cal G} = ({\cal V}, {\cal E})$\\ Number of cover nodes $K$}
\Output{Cover nodes $\{{\cal G}_i\}_1^K$}
Compute a graph clustering of ${\cal G}$ to $K$ clusters (graph cut or $K$-means on graph embeddings)\\
Assign each of the cluster to ${\cal G}_i,~i=1,\cdots,K$ \\
\While{condition not satisfied}{
Assign $H_j$ to ${\cal G}_i$ if $H_j$ is neighboring to ${\cal G}_i$ in ${\cal G}$ or within distance of $\epsilon$ to ${\cal G}_i$ in the embedding space\\
Build ${\cal K}$ from $\{{\cal G}_i\}_i^K$ \\
Check if ${\cal K}$ is connected \\
Compute homology group using \cite{Zomo2005} \\
Check if ${\cal H}_1({\cal K})$ is trivial \\
\If{Both conditions satisfied}{
Break
}
}
\end{algorithm}

To develop an algorithm based on Theorem~\ref{thm:consistency}, we proposed a greedy algorithm to construct ${\cal K}$ as in Algorithm \ref{alg:greedyK}. Note that a complex ${\cal K}$ is \emph{simply connected}, if  1) it is connected; 2) the 1-dimensional homology group ${\cal H}_1({\cal K})$ is trivial. Condition 1) could be easily verified by any graph traversal algorithm (BFS/DFS), and condition 2) can be verified computationally as in~\cite{Zomo2005}.



 
\section{Distributed Optimization}
\label{sec:formulation}

In this section, we introduce the proposed distributed formulation of recovering cycle-consistent maps from noisy pair-wise maps.

\subsection{Formulation}
Our formulation takes as input the pairwise base maps $\phi_{ij}$. We follow the state-of-the-art work on convex relaxation of enforcing cycle-consistent maps
\cite{Huang2013,Chen2014} to encode $\phi_{ij}$ into a data matrix $\overline{\mbx}_{ij}$. Following the common strategy for optimizing point-based maps, we relax $\phi_{ij}$ to be a partial map and/or soft map, i.e. $\overline{\mbx}_{ij}\in[0,1]^{m_i\times m_j}$, where $m_i$ denote the number of vertices in $H_i$.

Let $\overline{\mbx}_{{\cal V}}$ be the matching matrix that encodes pair-wise maps in its blocks, i.e.
\[
\overline{\mbx}_{\cal V} =
\begin{pmatrix}
\mathbf{I}  & \overline{\mbx}_{12} &  \cdots & \overline{\mbx}_{1n} \\
\overline{\mbx}_{21} & \mathbf{I} & \cdots & \vdots \\
\vdots & \cdots & \ddots & \vdots \\
\overline{\mbx}_{n1} & \cdots & \cdots & \mathbf{I} \\
\end{pmatrix}
\]

The goal of our formulation is to find a matrix $\mb{X}_{\cal V}$ that encodes cycle-consistent maps from the noisy input $\overline{\mbx}_{\cal V}$.
To achieve this goal, one observation is that the desired matching matrix $\mb{X}_{\cal V}$ is low-rank (c.f.~\cite{Huang2013}). Specififcally, we assume there is an universal object of size $m$, i.e. there are totally $m$ distinct entities for all the objects $H_i$'s in ${\cal V}$. For each object $H_i$, we have a latent map encoded by $\mb{A}_{H_i}\in \{0, 1\}^{m_i\times m}$, which maps a vertex from $H_i$ to an entity in the universal object. Let $\mb{A}_{\cal V}$ be a tall matrix that concatenates $\mb{A}_{H_i}$, i.e., $\mb{A}_{\cal V} = (\mb{A}_{H_1}^{T},\cdots, \mb{A}_{H_n}^{T})^{T}$. It is easy to see that the map matrix $\mb{X}_{\cal V}$ admits a low-rank factorization as $\mb{X}_{\cal V} = \mathbf{A}_{\cal V} \mathbf{A}_{\cal V}^\top$. In~\cite{Huang2013, Zhou2015}, the authors use this property to develop robust algorithms for recovering $\mb{X}_{\cal V}$ from noisy input maps.

Without losing generality. we assume $h(\mb{X}_{\cal V})$ is an objective function that measures the quality of a set of cycle-consistent maps encoded by $\mb{X}_{\cal V}$, e.g., it promotes the low-rankness of $\mb{X}_{\cal V}$. The precise expression of $h(\mb{X}_{\cal V})$ will be discussed later. Our distributed formulation is given by
\begin{equation}\label{eqn:dsdp}
\begin{array}{ll}
\min & \sum_i h\left(\mb{X}_{{\cal V}_i}\right) \\
\textnormal{s.t.} & \mb{X}_{{\cal V}_{i\cap j}^i} = \mb{X}_{{\cal V}_{i\cap j}^j}, \forall (i,j) \in {\cal E},
\end{array}
\end{equation}
where $\mb{X}_{{\cal V}_{i\cap j}^i}$ is the matching matrix of ${\cal V}_{i\cap j}$ in ${\cal V}_i$, i.e. a sub-matrix of $\mb{X}_{{\cal V}_i}$ by picking blocks that belong to the matching graphs in ${\cal V}_{i\cap j}$. Each $h(\mb{X}_{{\cal V}_i})$ indicates local consistency in ${\cal V}_i$, and the condition that $\mb{X}_{{\cal V}_{i\cap j}^i} = \mb{X}_{{\cal V}_{i\cap j}^j}$ will guarantee that the overlapping subgraph are consistent. In such a manner, the consistency condition will propagate through the overlapping sub-graph to each component ${\cal V}_i$ conceptually similar to our proof of Theorem \ref{thm:consistency}.

In the state-of-the-art methods of~\cite{Huang2013} and~\cite{Zhou2015}, the authors proposed different formulations of objective function $h(\mb{X}_{\cal V})$. We will use the formulation described in~\cite{Zhou2015} to demonstrate our framework, because of its competent performance and superior computational speed.


As in \cite{Zhou2015}, $h(\mb{X}_{\cal V})$ can be written as
\begin{equation}\label{eqn:gsdp}
\begin{array}{ll}
\min & \langle \mb{W}_{\cal V}, \mb{X}_{\cal V}\rangle + \lambda\|\mb{X}_{\cal V}\|_* \\
\textnormal{s.t.} & \mb{X}_{\cal V} \succeq 0, \\
& \mb{X}_{{\cal V}(ii)} = \mathbf{I}_{m_i}, \forall i\\
& \mb{X}_{{\cal V}(ij)} = \mb{X}_{{\cal V}(ji)}^\top, \forall i\neq j\\
& 0 \leq \mb{X} \leq 1\\
\end{array}
\end{equation}
where $\langle\cdot,\cdot\rangle$ is the matrix inner product, \textcolor{black}{$\|\cdot\|_*$ is the matrix nuclear norm}, and \textcolor{black}{$\mb{W}_{\cal V} =  \alpha\mathbf{1} - \overline{\mbx}_{\cal V}$}, and $\mathbf{1}$ denote the matrix whose elements are $1$. The purpose of adding constant $\alpha$ is to impose a $L_1$ constraint on \textcolor{black}{$\mb{X}_{\cal V}$} to promote sparsity.
\textcolor{black}{We use $\mb{X}_{{\cal V}(ij)}$ to denote the $(i,j)^\textnormal{th}$ block of the block matrix $\mb{X}_{\cal V}$}.
As has been shown in \cite{Zhou2015}, the constraint $\mb{X}_{\cal V} \succeq 0$ may be relaxed for a sufficiently large $\lambda$. Let ${\cal C}_i$ encode the convex set induced by the constraints for ${\cal V}_i$, we could then simplify the formulation of our distributed problem as
\begin{equation}\label{eqn:dsdp2}
\begin{array}{ll}
\min &\textcolor{black}{\sum_i\left(\langle\mb{W}_{{\cal V}_i}, \mb{X}_{{\cal V}_i}\rangle + \lambda\|\mb{X}_{{\cal V}_i}\|_*\right)} \\
\textnormal{s.t.} & \mb{X}_{{\cal V}_i} \in {\cal C}_i\\
& \mb{X}_{{\cal V}_{i\cap j}^i} = \mb{X}_{{\cal V}_{i\cap j}^j}, \forall (i,j) \in {\cal E}
\end{array}
\end{equation}


\section{Alternating minimization}
\label{sec:solver}
\subsection{Algorithms}
\label{ssect:alg}
The nuclear norm minimization in \eqref{eqn:dsdp2} can be efficiently optimized using recent results on low-rank matrix recovery, which directly enforce low-rank decompositions $\mb{X}_{{\cal V}_i} = \mb{A}_{{\cal V}_i}\mb{B}_{{\cal V}_i}^\top$ (c.f.~\cite{Cabral2013, Hastie2015, Zhou2015}). Here $\mb{A}_{{\cal V}_i}$ and $\mb{B}_{{\cal V}_i}$ are latent variables. According to~\cite{Recht2010}, we can write the nuclear norm as
\[
\|\mb{X}\|_* = \min_{\mb{A}, \mb{B}: \mb{A}\mb{B}^\top = \mb{X}} \frac{1}{2}\left(\|\mb{A}\|_F^2 + \|\mb{B}\|_F^2\right).
\]
To make the notations uncluttered, we will shorten $\mb{X}_{{\cal V}_i}$ and $\mb{X}_{{\cal V}_{i\cap j}}^i$ as $\mb{X}_i$ and $\mb{X}_{ij}$, respectively.
Moreover, let $\mb{E}_{ij}$ denote the selection matrix to extract the part of $\mb{X}_i$ that belongs to the set of ${\cal V}_i\cap{\cal V}_j$, i.e.
$\mb{X}_{ij} = \mb{E}_{ij}^\top\mb{X}_i\mb{E}_{ij}$. With this setup, the condition on intersection consistency
becomes $\mb{E}_{ij}^\top\mb{X}_i\mb{E}_{ij} = \mb{E}_{ji}^\top\mb{X}_j\mb{E}_{ji}$.

We then finalize our formation of the problem in \eqref{eqn:dsdp} as
\begin{equation}\label{eqn:dsdp-final}
\begin{array}{ll}
\min & \sum_i\left(\langle \mb{W}_i, \mb{X}_i\rangle + \frac{\lambda}{2}\|\mb{A}_i\|_F^2 + \frac{\lambda}{2}\|\mb{B}_i\|_F^2\right) \\
\textnormal{s.t.} & \mb{X}_i = \mb{A}_i\mb{B}_i^\top, \\
& \mb{E}_{ij}^\top\mb{X}_i\mb{E}_{ij} = \mb{E}_{ji}^\top\mb{X}_j\mb{E}_{ji}, \\
& \mb{X}_i \in {\cal C}_i,
\end{array}
\end{equation}
In all our experiments, we set $\alpha=0.1$, $\lambda = 50$, $\mu = 64$, and $\beta=1$.

We apply ADMM to solve \eqref{eqn:dsdp-final}. The solver
is summarized in Algorithm \ref{alg:dsdp}. In particular, $\mb{Y}_i$ and $\mb{Z}_{ij}$ are dual variables. The constraints on $\mb{X}$ are handled implicitly and updated in the alternating algorithm. The ADMM algorithm updates primal variables by minimizing ${\cal L}$ and then applies gradient descent to update the dual variables.
\begin{algorithm}[ht]
\caption{Distributed Graph Matching via ADMM}
\label{alg:dsdp}
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{Initial pairwise maps $\overline{\mbx}_i$}
\Output{Consistent matches $\mb{X}_i$}
Initialize $\mb{A}_i$, $\mb{B}_i$ randomly, and set $\mb{Y}_i$, $\mb{Z}_{ij}$ to be $\mathbf{0}$ \\
$\mb{W}_i = \alpha\mathbf{1} - \overline{\mbx}_i$ \\
\While{not converged}{
\tcc{inner-node update}
$\mb{A}_i \leftarrow (\mb{X}_i + \frac{1}{\mu}\mb{Y}_i)\mb{B}_i(\mb{B}_i^\top\mb{B}_i + \frac{\lambda}{\mu}\mathbf{I})^\dagger$ \\
$\mb{B}_i \leftarrow (\mb{X}_i + \frac{1}{\mu}\mb{Y}_i)\mb{A}_i(\mb{A}_i^\top\mb{A}_i + \frac{\lambda}{\mu}\mathbf{I})^\dagger$ \\
$\mb{X}_i \leftarrow {\cal P}_{{\cal C}_i}\left(\mb{X}_{i0}\right)$ \\
$\mb{Y}_i \leftarrow \mb{Y}_i^k + \mu\left(\mb{X}_i - \mb{A}_i\mb{B}_i^\top\right)$ \\
$\mb{Z}_{ij} \leftarrow \mb{Z}_{ij}^k + \beta\left(\mb{M}_{i\rightarrow j}^{k+1} - \mb{M}_{j\rightarrow i}^{k+1}\right)$ \\
\tcc{inter-node information exchange}
node $j$ send $\mb{M}_{j\rightarrow i} = \mb{E}_{ji}^\top\mb{X}_j\mb{E}_{ji}$ to node i
}
Round $\mb{X}_i$ with a threshold of 0.5.
\end{algorithm}
Moreover, $\mb{A}_i$ and $\mb{B}_i$ admit closed-form solution via solving least-squares. Moreover, $\mb{X}_{i0}$ is the solution to the linear equation
\begin{multline}
\mu\mb{X}_i + 2\beta\sum_j\mb{E}_{ij}\mb{E}_{ij}^\top\mb{X}_i\mb{E}_{ij}\mb{E}_{ij}^\top = \mu\mb{A}_i\mb{B}_i^\top - (\mb{W}_i + \mb{Y}_i) \\
+ \sum_j\mb{E}_{ij}(2\beta\mb{M}_{j\rightarrow i}^k - \mb{Z}_{ij} + \mb{Z}_{ji})\mb{E}_{ij}^\top. \nonumber
\end{multline}
Furthermore, the update on $\mb{X}_i$ requires a projection onto the convex set ${\cal C}$, ${\cal P}_{\cal C}(\cdot)$, i.e. ${\cal P}_{\cal C}(\mb{X}_0)$ is the solution to the problem
\[
\min_{\mb{X}\in{\cal C}}\|\mb{X} - \mb{X}_0\|_F^2.
\]
This is essentially a linear programming problem, and can be solved efficiently. We refer to the supplementary material for a derivation.

The key point in Algorithm \ref{alg:dsdp} is the separation of inner-node update and inter-node information exchange. It can be seen that all the matrix computations are done in each node separately, which indicates a distributed computation. While after each iteration, adjacent nodes will need to exchange information by passing messages. Namely, for node ${\cal V}_j$ and all its neighboring nodes ${\cal V}_i$, a message will be sent in the form of $\mb{M}_{j\rightarrow i} = \mb{E}_{ji}^\top\mb{X}_j\mb{E}_{ji}$. Note that there is no overhead on generating these messages. Recall that $\mb{E}_{ji}$ is just a sub-block extraction matrix, and the way $\mb{M}_{j\rightarrow i}$ is computed is simply by extracting the sub-block of the matching matrix $\mb{X}_j$ that belongs to the intersection ${\cal V}_i\cap{\cal V}_j$. Since all the computations are indeed done on each node individually, the proposed algorithm is essentially completely distributed, and the only add-on is a syncing stage.


\subsection{Complexity}
The computational complexity of Algorithm \ref{alg:dsdp} is dominated by matrix multiplication. In our approach, the per-iteration complexity is controlled by the leading node in ${\cal G}$, i.e. $O\left(\max_i(\sum_{H_j\in{\cal V}_i}m_j)^2\max_i(|{\cal V}_i|)\right)$, where $|{\cal V}_i|$ is the total number of distinct entities among all the objects in ${\cal V}_i$. In contrast, the per-iteration complexity of~\cite{Chen2014} and~\cite{Zhou2015} is $O((\sum_im_i)^3)$ and $O((\sum_im_i)^2m)$ respectively. Furthermore, in our experiments, we found the total number of iterations to converge for our method is comparable to that of ~\cite{Zhou2015}.


\section{Experiments}
\label{sec:exp}
\subsection{Simulation}

In this section, we perform experimental evaluation using synthetic datasets.

We followed the same experimental setup as in \cite{Chen2014, Zhou2015}. Given an optimized matching matrix $\mb{X}^*$ and the ground truth mapping $\mb{X}_g$, we access the quality of $\mb{X}^*$ by measuring the intersection over union (or IOU) score:
\[
1 - \frac{|\tau(\mb{X}^*)\cap\tau(\mb{X}_g)|}{|\tau(\mb{X}^*)\cup\tau(\mb{X}_g)|}
\]
where $\tau(\cdot)$ denotes the mapping induced from the matching matrix, and $|\cdot|$ denotes the set size. Note in our distributed setting, we could only partially recover $\mb{X}^*$ given $\mb{X}_i^*$. Therefore, our ground truth setting is also different from \cite{Chen2014, Zhou2015} in this regard.

\subsubsection{Matching Errors}
\label{sssect:merror}

In the first experiment, we aimed to evaluate the matching performance between our algorithm, DMatch, and the global algorithm, MatchALS, as in \cite{Zhou2015}. We considered the following model to generate the testing examples. The total number of graphs is denoted by $n$. The size of the universe is fixed at $r=20$ points. In each graph, a point is randomly observed with a probability $\rho_0$. We simulated error corruption by randomly removing true mapping and adding false ones with a corruption rate $\rho_e$.

We considered two settings in our experiments. In the first setting, we constructed our cover graph by making a sparse three way tree. This was done by randomly selecting a subset of ${\cal V}$ as a common intersection ${\cal V}_c$ and then split the rest evenly into the three cover nodes ${\cal V}'_1,{\cal V}'_2,{\cal V}'_3$. As a consequence, each cover node is equal to ${\cal V}_i = {\cal V}_c\cup{\cal V}'_i$.
In the second setting, we increased the overlap density by circularly adding one more split to each cover, i.e. ${\cal V}_i = {\cal V}_c\cup{\cal V}'_i\cup{\cal V}'_{i+1}$. We compared DMatch to MatchALS by varying the parameters $\rho_0$, $\rho_e$, and $n$. For both algorithms, we used $m=2r$ and $\lambda=50$.

\begin{figure}[ht]
\centering
\begin{subfigure}{.96\linewidth}
\begin{subfigure}{.32\linewidth} \centering \includegraphics[width=0.96\linewidth]{sparse_match_error_n}
\vspace{-10pt}\center{\footnotesize DMatch-sparse}
\end{subfigure}
\begin{subfigure}{.32\linewidth} \centering \includegraphics[width=0.96\linewidth]{dense_match_error_n}
\vspace{-10pt}\center{\footnotesize DMatch-dense}
\end{subfigure}
\begin{subfigure}{.32\linewidth} \centering \includegraphics[width=0.96\linewidth]{m_match_error_n
\vspace{-10pt}\center{\footnotesize MatchALS}
\end{subfigure}
\caption{number of objects $n$ + input error rate $\rho_e$}
\end{subfigure}
\centering
\begin{subfigure}{.96\linewidth}
\begin{subfigure}{.32\linewidth} \centering \includegraphics[width=0.96\linewidth]{sparse_match_error_obs}
\vspace{-10pt}\center{\footnotesize DMatch-sparse}
\end{subfigure}
\begin{subfigure}{.32\linewidth} \centering \includegraphics[width=0.96\linewidth]{dense_match_error_obs}
\vspace{-10pt}\center{\footnotesize DMatch-dense}
\end{subfigure}
\begin{subfigure}
{.32\linewidth} \centering \includegraphics[width=0.96\linewidth]{m_match_error_obs}
\vspace{-10pt}\center{\footnotesize MatchALS}
\end{subfigure}
\caption{observation ratio $\rho_0$ + input error rate $\rho_e$}
\end{subfigure}
\caption{Matching error comparison. Darker color means lower matching error. The fixed parameter is set to $\rho_0=0.6$ and $n = 50$ respectively.}
\label{fig:merror}
\end{figure}

Figure~\ref{fig:merror} shows matching errors under various configurations, for both DMatch and MatchALS. In general, lowering input error and increasing observation ratio or increasing the total number of objects will improve the matching performance, i.e. with a lower matching error. In addition, we can see that increasing the overlap between cover nodes would have a positive impact on the recovery (comparison between Figure \ref{fig:merror}(a) and Figure \ref{fig:merror}(b)).

Furthermore, in a comparison between Figure \ref{fig:merror}(b) and Figure \ref{fig:merror}(c), we can see that when the cover is dense enough, i.e. the size of the overlaps are sufficiently large, the matching error would approach that of MatchALS, which is the global recovery.

\subsubsection{Graph Covers}
In the second experiment, we aimed to understand more on the effects of graph covers. Specifically, we construct a ground truth graph cover by selecting a sparse cover as in Section \ref{sssect:merror}. For every pair of graphs within the same cover node, we set the error rate to be $\rho_{in}$, and for every pair between different cover node, we set the error rate to be $\rho_{out}$. The experiment is then conducted by comparing DMatch to 1) using the ground truth cover and 2) using a randomly constructed cover.

\begin{figure}[ht]
\centering
\begin{subfigure}{.49\linewidth}
\centering
\begin{subfigure}{.49\linewidth} \centering \includegraphics[width=0.96\linewidth]{sparse_gt_cover_merror}
\centering{\footnotesize sparse}
\end{subfigure}
\begin{subfigure}{.49\linewidth} \centering \includegraphics[width=0.96\linewidth]{dense_gt_cover_merror}
\centering{\footnotesize dense}
\end{subfigure}
\caption{\footnotesize ground truth cover}
\end{subfigure}
\centering
\begin{subfigure}{.49\linewidth}
\centering
\begin{subfigure}{.49\linewidth} \centering \includegraphics[width=0.96\linewidth]{sparse_rand_cover_merror}
\centering{\footnotesize sparse}
\end{subfigure}
\begin{subfigure}{.49\linewidth} \centering \includegraphics[width=0.96\linewidth]{dense_rand_cover_merror}
\centering{\footnotesize dense}
\end{subfigure}
\caption{\footnotesize randomly constructed cover}
\end{subfigure}
\caption{Effects on recovery from the construction of graph cover. Darker color means lower matching error. The fixed parameter is set to $\rho_0=0.6$ and $n = 50$.}
\label{fig:gcover_error}
\end{figure}

Figure~\ref{fig:gcover_error} shows the experimental results. In Figure~\ref{fig:gcover_error}(a), we can see that using the ground truth cover, the matching performance does not depend much on $\rho_{out}$, because we explicitly disregarded any information from $\rho_{out}$. While on the randomly constructed cover, the two error rates are mixed up. Specifically, the results favor more on small $\rho_{in}$ and small $\rho_{out}$ at the same time. This situation, however, starts to change when the cover becomes denser, i.e. the increasing overlap between cover nodes. The dependency on $\rho_{out}$ disappears as shown in Figure \ref{fig:gcover_error}(b). One potential explanation is that with the increasing overlap between cover nodes, the portion of out-node pairs become smaller, as well as the resultant portion of induced error from them. As a consequence, the mixed error rate is dominated by $\rho_{in}$. In addition, we could also see that with denser cover, the algorithm is more error tolerant. This comes with a trade-off that on average the size of each cover becomes bigger and the computational cost also increases.

\subsection{Real Experiments}

\subsubsection{CMU House Sequence}
In this part of the experiment, we intend to test the scalability of our distributed algorithm. We use the CMU House sequence\footnote{\url{http://vasc.ri.cmu.edu//idb/html/motion/house/index.html}} as the testing images. This sequence has been widely used to test different graph matching algorithms. It consists of 110 frames, and there are 30 feature points labeled consistently across all frames. We used the geometry based constraint in pairwise matching as is done in \cite{Hu2014}. To construct a valid cover complex ${\cal K}$, we first computed all pairwise matches and for each match, the result is encoded using a binary matrix. We then built a fully connected matching quality graph, where each vertex represents an image, and the edge weight represents the matching score associated with each image pair. Since we knew the sequence was roughly generated by moving a camera in a single dimension, we computed the Laplacian embedding of the vertices onto 1 dimensional space using Fielder vector, and then we applied Algorithm \ref{alg:greedyK} to build the cover complex.
We compared our DMatch algorithm with the global methods, MatchALS~\cite{Zhou2015}, MatchLift~\cite{Chen2014}, Spectral ~\cite{Pachauri2013}, and the distributed method described in~\cite{Leonardox2017}. Besides running global algorithms on fully connected map graph, we also ran them on the sparse map graph induced by our cover complex, marked with different $K$ values. The algorithm was implemented in a single laptop with Intel Core i7 2.0GHz CPU and 16GB DDR3 Memory. We measured the time used in each cover node separately and then took the max as the total computational time per iteration, where we assumed the cost for messages passed between adjacent cover nodes was negligible.

\begin{table}[ht]
\centering
  \begin{tabular}{ r r r r}
    \hline \hline
      & Error Rate & Iterations & Time \\ \hline
    Original & 0.1445 & - & - \\
    MatchALS & 0.1031 & 266 & 98.8 \\
    MatchLift & 0.1027 & 1000 & 3791.1 \\
    Spectral & 0.1277 & - & \bf{0.6} \\
    MatchALS ($K=4$) & 0.0161 & 380 &103.5 \\
    MatchALS ($K=6$) & 0.0648 & 1000 & 268.2 \\
    MatchLift ($K=4$) & \bf{0} & 1000 &4066.1 \\
    MatchLift ($K=6$) & \bf{0} & 1000 & 3972.3 \\
    DMatch ($K=4$) & \bf{0} & 203 & 28.9 \\
    DMatch ($K=6$) & \bf{0} & \bf{150} & 7.6 \\
    Consensus ($K=4$) & \bf{0} & 150 & 28.3 \\
    Consensus ($K=6$) & 0.0071 & 150 & 18.7 \\
    \hline \hline
  \end{tabular}
\caption{The error rate and the total computational time (seconds) on CMU House sequence.}
\label{tab:scalability}
\end{table}

Table~\ref{tab:scalability} shows the matching accuracy, timing and iterations used in these algorithms. In our experiments, Spectral method appeared to be the fastest method among all the methods. On the other hand, it has the largest reported error. The consensus algorithm is not error-driven, hence we used a preset 150 iterations to match the lowest number of iterations among all experiments. In general, distributed algorithms used less iterations to converge, and achieved at least an order of magnitude speed-up compared with global methods, while maintaining an error rate of {\bf 0} (almost 0 for Consensus algorithm with $K=6$). It can also be seen that for distributed algorithms once we increased the number of clusters from $K=4$ to $K=6$, both the number of iterations and total computational time decreased. This, on the other hand, proved our complexity analysis, since increasing the number of clusters would in general reduce the number of vertices in each cluster. Although both DMatch and Consensus algorithms achieved similar results in terms of accuracy and time, the latter requires knowledge of the number of universal entities and has limitations dealing with partial matches. Another interesting observation is after convergence, the error of MatchALS on the sparse map graph induced by our cover complex is reduced, when compared with that on the fully connected map graph. One explanation is that the graph cover structure grouped together images that have high pairwise matching quality and explicitly disregard any pairwise matches that are of low quality (covers are joint normal), and as a consequence, it is robust against noisy pair-wise matches. In addition, we also extended the experiments for $K=8,10$ and found the running time was not reduced significantly, in comparison with that from $K=4$ to $K=6$. Since we need to have a reasonable amount of overlaps between clusters to pass matching information around, increasing the number of clusters does not necessarily reduce the maximum size of the clusters.

\begin{figure*}[t]
\begin{subfigure}{.32\linewidth} \centering \includegraphics[width=0.96\linewidth]
{graff_curve}
\caption{Graffiti}
\end{subfigure}
\begin{subfigure}{.32\linewidth} \centering \includegraphics[width=0.96\linewidth]
{bikes_curve}
\caption{Bikes}
\end{subfigure}
\begin{subfigure}{.32\linewidth} \centering \includegraphics[width=0.96\linewidth]
{leuven_curve}
\caption{Leuven}
\end{subfigure}
\caption{The performance curve on Graffiti, Bikes and Leuven datasets. The $y$-axis is the correct match ration and the $x$-axis is the threshold value over the image width. We compare DMatch (red solid) with MatchALS \cite{Zhou2015} (blue dashed), MatchLift \cite{Huang2013}, Spectral \cite{Pachauri2013}, and original pairwise matching (black dotted).}
\label{fig:error_curve}
\end{figure*}

\begin{figure*}[t]
\begin{subfigure}{.3\linewidth} \centering \includegraphics[height=6cm]{graffiti_matches}
\caption{Graffiti}
\end{subfigure}
\begin{subfigure}{.34\linewidth} \centering \includegraphics[height=6cm]{bikes_matches}
\caption{Bikes}
\end{subfigure}
\begin{subfigure}{.34\linewidth} \centering \includegraphics[height=6cm]{light_matches}
\caption{Leuven}
\end{subfigure}
\caption{Example of matching results. The bottom one is from DMatch, the middle one is from MatchALS, and the top one is from original pairwise matching respectively. Yellow lines encode the correct matches, while blue lines are for wrong ones.}
\label{fig:matchings}
\end{figure*}

\subsubsection{Graffiti datasets}
In this experiment, we followed the procedure described in \cite{Zhou2015}. We used the benchmark datasets from Graffiti datasets \footnote{\url{http://www.robots.ox.ac.uk/~vgg/data/data-aff.html}}. In each dataset, there are 6 images of a scene with various image transformations, including viewpoint change, blurring, and illumination variation etc.

To construct an affinity score matrix $\overline{\mbx}$, we employed the same procedure as in \cite{Zhou2015} for comparison purpose. We first detected 1000 affine covariant feature \cite{Miko2005} points in each image of the dataset and computed their SIFT \cite{lowe2004} descriptors using VLFeat library \cite{vlfeat}. The affinity scores were computed as the inner product between every pair of detected feature points on each pair of images.
We excluded apparent mismatches by keeping only affinity scores that are above the threshold $0.7$. Furthermore, any potential matches that are indistinguishable was removed, i.e. if the first and the second top matches were below the ratio threshold $1.1$, the candidate point was removed. Finally, any feature point that has only one candidate match
in the dataset was also excluded.

In a comparison, to construct our cover graph,
we first built a matching quality graph, using the matching score as the edge weight and
used the Fiedler vector of the graph laplacian as the embedding and applied Algorithm \ref{alg:greedyK} to build the cover complex.

To evaluate the performance, we used the ground truth homography matrix given in the dataset, and adopted the procedure used in \cite{Chen2014}. For a testing point, we calculated the true correspondence using homography and compared with the matched correspondence. If they were within a predefined distance threshold, we deemed the matching is correct, and otherwise, wrong. Then we swept along the threshold dimension to draw an error curve that is dependent on the threshold chosen.

We tested our DMatch algorithm against MatchALS~\cite{Zhou2015}, MatchLift~\cite{Chen2014}, Spectral~\cite{Pachauri2013} and the original pairwise matchings. We ignored the Consensus algorithm \cite{Leonardox2017} as it cannot explicitly handle partial matches.
Figure~\ref{fig:error_curve} shows the curve for three datasets, Graffiti, Bikes, and Leuven. Note that DMatch will not give a full pairwise matching between images, instead, we only have a matching when the two images belong to the same cover node. Therefore, we computed the one-hop composite match between image pairs across different cover nodes\footnote{For each image $i$ and $j$ not in the same cover node, we loop through all $k\neq i,j$ and accumulate the composite matchings from $i\rightarrow k$ and $k\rightarrow j$.}. From the performance curve, we can see that our DMatch performs very similar to the best global methods in all datasets as shown in \cite{Zhou2015}. In another word, DMatch achieved performance gains without loss of matching quality.

In Figure \ref{fig:matchings}, we show the example matches between the first and the fourth image for each dataset. The bottom match is DMatch, the middle is MatchALS and the top one is the original pairwise map. Clearly, our matching shows at least as good as the results of MatchALS, where both corrected mismatches (reduced blue lines) and increased correct matches (denser yellow lines).

In our implementation
we notice that the total number of iterations to converge for both DMatch and MatchALS are roughly the same (around 60 iterations).

\section{Conclusion}
In this paper, we introduced a scalable framework for establishing consistent matches across multiple graphs in a distributed manner.
We showed how to use our framework to extend state-of-the-art global methods. By running an iterative optimization algorithm locally and exchange information in every iteration, our framework would achieve local and global consistent matching at the same time. Furthermore, we theoretically proved the sufficient conditions under which locally consistent matching would guarantee global consistency. In our experiments, we showed that in practice, the assumptions and the conditions in our theorem could be relaxed without sacrificing performance. In addition, our proposed distributed framework achieved order of magnitude improvements in speed. We believe this is a very important first step for large scale exploration of images for object matching as well as building 3D object models from crowd-sourced collections. Future work includes matching large collection of different deformable objects that have high similarity and enough variance, e.g. a collection of different dogs or cats.

\vspace{10pt}
\noindent{\bf Acknowledgements:} The authors wish to thank the support of NSF grants DMS-1546206 and DMS-1700234, a Google focused research award and a gift from Adobe Research.




{\small
\bibliographystyle{ieee}
