\section{Introduction}
The \emph{assignment problem} finds an assignment, or matching, between two finite sets
$U$ and $V$, each of cardinality $n$, such that the total cost of all matched pairs is minimized. The assignment problem can also be generalized to finding matchings between more than two sets. 
This is a fundamental problem in computer science and has been motivated by a wide gamut of research areas spanning diverse areas such as structural biology~\cite{singer2011three}, protein structure comparisons in bioinformatics~\cite{zaslavskiy2009global}, and computer vision~\cite{conte2004thirty}. Computer vision especially boasts a broad range of applications that  include object matching, image registration~\cite{shen2002hammer}, stereo matching~\cite{goesele2007multi}, shape matching~\cite{petterson2009exponential,berg2005shape}, structure from motion (SfM)~\cite{szeliski2010computer}, and object detection~\cite{jiang2011linear}, to name a few.

Various assignment approaches can 
broadly be classified as those that find a bijective assignment in the form of a permutation matrix by posing the problem as a \emph{linear assignment problem (LP)} versus ones that solve a \emph{quadratic assignment problem (QAP)} via \emph{graph matching}, where each graph's nodes represent the objects and the edges encode their corresponding distances; the goal of QAP then is to find node-wise correspondences between the graphs so that the overall discrepancy between their corresponding edge-wise counterparts is minimized and the overall \emph{relational structure} is best preserved.

\emph{Partial assignment} implies that only subsets of $U$ and $V$ can actually be assigned to each other successfully. 
This phenomenon is of particular interest to applications where either objects are absent due to incomplete observations, undergo deformations, and/or the objects in question cannot clearly be disambiguated because the objects in question along with their related objects are embedded in clutter. This variant of the assignment problem is widely accepted as a formidable challenge.

\begin{figure}[t!]
	\centering
	\includegraphics[width=\columnwidth]{Experiments/drawing.png}
	\caption{Matching cliques of two houses with $1$-cliques (red), $2$-cliques (green), $3$-cliques (blue). Violet and yellow lines show the matchings of $d$-cliques where $d=2,3$ respectively.}
	\label{fig:drawing}
\end{figure}

Although graph matching methods were found to be instrumental, they too perform poorly when faced with non-similar geometric transformations or transformations that produce degenerate triangulations. This is attributed to assigning weights to only node and edge assignments, while ignoring the interplay of higher-order connections/relations. For example, using triplet weights can alleviate the above mentioned problem to a very large extent by defining a measure invariant to scale and other transformations~\cite{Chertok2010}.

Motivated by the aforementioned observations and inspired by Kahle~\cite{kahle2006}'s work on combinatorial topological models like the \emph{random clique complex}, we focus our attention 
to matching higher-order components between two sets of points in the setting of some points \emph{missing completely at random}. We pose our assignment problem as finding a matching between two sets of points, each represented as a random clique complex, which is a higher-order analogue of random graphs. Figure~\ref{fig:drawing} illustrates such a matching of cliques of corresponding dimensionality, between two different scenes (taken from different camera angles) of the same house. 
Given an Erdős-Rényi (ER) graph, its \emph{clique complex} is the simplicial complex with all complete subgraphs (i.e., cliques) as its \emph{faces}. The Erdős-Rényi graph forms the $1$-skeleton of the random clique complex, where the cliques have at most a dimension of $2$, i.e., edges in the graph.
The clique topology of a random adjacency matrix is analogous to its eigenvalue spectrum, as it provides a set of invariants that help detect structure~\cite{giusti2015}. This probabilistic and combinatorial framework of random clique complexes allows us to further study the assignment problem under various assumptions of the underlying distribution of the matrix entry distributions, 
its robustness to missing values, and its asymptotic behavior for large-scale cases.

\textbf{Contributions:} We present the following contributions.
\begin{enumerate}
	\item To the best of our knowledge, our proposed approach is a first attempt to formulate 
	\emph{higher-order matching} between two sets of points, given partial or incomplete information, as a matching between two \emph{random clique complexes}. We also propose an efficient matching algorithm and study both its time and storage complexity.
	
	\item (i) We provide new bounds on the concentration inequality of eigenvalues of the QAP trace formulation for random symmetric matrices, (ii) we give tighter concentration inequality bounds on the largest eigenvalue for the Lawler QAP formulation on random matrices, in the context of \emph{affinity matrices} that are used by some earlier works.  Furthermore, we theoretically analyze and discuss the robustness of affinity-matrix based schemes to missing points, and
	(iii) we perform asymptotic analysis on the worst to best case ratio of a QAP solution for our higher-dimensional clique adjacency matrices in the \emph{clique percolation regime}~\cite{Bollobas2009}, 
	where the entries follow a \emph{Poisson} distribution.
	  
	\item Finally, we present a comprehensive empirical study that compares our method's matching accuracy to that of a diverse set of matching approaches~\cite{ZhouD16,zhou2013deformable,cho2010reweighted,feizi2016spectral,Leordeanu2005,cour2007balanced,pachauri2013solving,gold1996graduated,Kuhn1955,leordeanu2009integer,zass2008probabilistic,li2013object,duchenne2011tensor}. 
	We conducted our experiments on both synthetic and well-known hard real-world datasets that span across affine/non-affine transformations, severe occlusions, and clutter. Our study reveals much better accuracy for the popular datasets against several of the state-of-the-art matching methods.
\end{enumerate}

\section{Matching Random Clique Complexes}
\label{sec:match_rcc}
We consider the problem of capturing higher-order feature groups among landmark points in an image by representing them as a \emph{random clique complex} (RCC) and then using these RCCs to match two sets of groupings from two different images. We begin this section by describing the construction of a random clique complex, followed by our proposed method of matching two RCCs, and we finally analyze the runtime and storage complexity of our algorithm. 

\subsection{Structure of a Random Clique Complex}
\label{subsec:rcc}
We begin with general definitions pertaining to the structure of simplicial complexes and then accordingly adapt these definitions to our domain of random graphs to build random clique complexes.

Let $G(n,p)$ be an Erdős-Rényi graph with a set of $n$ vertices denoted by $V$, whose edges $\{v,v'\} \in {V \choose 2}$, are i.i.d Bernoulli($p$) distributed. Recall, that a $k$-clique in $G(n,p)$
is a \emph{complete subgraph} that comprises of $k$ vertices and ${k \choose 2}$ edges.
Here onwards, for ease of notation, we will denote $G(n,p)$ as $G$. Given any affinely independent set $V=\{v_i\}_{i=0}^k$ of $(k+1)$ points in $\mathds{R}^n$, the $k$-simplex $\sigma^{(k)}$ is the \emph{convex hull} of $V$,
i.e., it is the set of all points of the form $w_0 v_0 + \dots + w_k v_k$, where 
$\sum_{i=0}^{k} w_i = 1$ and $w_i \geq 0$ for all $i$. 
If we imagine the vertices of $G$ embedded generically in $\mathds{R}^n$, then each $(k+1)$-clique consisting of $k+1$ vertices is represented by a $k$-dimensional simplex $\sigma^{(k)}$ in our random clique complex. For example, a $2$-clique (edge) and a $3$-clique (triangle) in $G$ is represented as $\sigma^{(1)}$ and $\sigma^{(2)}$, respectively.

Given $0\leq i \leq k$, the $i$-th face $f_i$ of $\sigma^{(k)}$ is the subspace of points that satisfy $w_i=0$; it is the $(k-1)$-simplex $\sigma^{(k-1)}$ whose vertices are all those of $\sigma^{(k)}$, except the $i$-th vertex. In other words, when 
$\sigma^{(k)}$ is a clique of $G$, then all its subsets are also cliques and hence 
considered faces of $\sigma^{(k)}$. For example, a $3$-clique (triangle) has three $2$-cliques (edges) in it.

With the aforementioned definitions in mind, we define our random clique complex 
$\mathcal{X}(G)$ as the set of all cliques in $G$ such that
$\mathcal{X}(G) = \{ \sigma \in [n] \mid \sigma \text{ is a clique of }G  \}$. We denote a set of $(k+1)$-cliques as $\mathcal{X}_k(G)$. Additionally, $\mathcal{X}(G)$ also satisfies the following conditions of a simplicial complex: (i) Any face in $\mathcal{X}(G)$ is also a simplex in $\mathcal{X}(G)$ and (ii) the intersection of any two simplexes $\sigma_i, \sigma_j$ is a 
face (lower dimensional clique) of both $\sigma_i$ and $\sigma_j$.

The faces of $\sigma^{(k)}$ are \emph{copies} of $\sigma^{(j)}$ for $j<k$, which are glued together inductively. The $k$-skeleton of $\mathcal{X}(G)$, for $k \in \mathds{N}$, is defined as the following quotient space
\begin{align*}
\mathcal{X}^{(k)}(G) :=  \left( \mathcal{X}^{(k-1)}(G) \cup \coprod_{\sigma: \textit{dim } \sigma=k}  \sigma^{(k)}  \right) \Bigg/ \sim
\end{align*}

where $\sim$ is the equivalence relation that identifies faces of $\sigma^{(k)}$
to the corresponding faces of $\sigma \in \mathcal{X}^{(j)}(G)$ where $j<k$. 
Finally, $\mathcal{X}(G)=\cup_{k=0}^\infty \mathcal{X}^{(k)}(G)$.

\textbf{$k$-skeleton as adjacency matrix}: Given a random graph $G$ and its $k$-skeleton 
$\mathcal{X}^{(k)}(G)$ that contains all its $(k+1)$-cliques, we follow the idea from Bollob{\'{a}}s et. al.~\cite{Bollobas2009}, to 
represent $\mathcal{X}^{(k)}(G)$ as an adjacency matrix $G^{(k,l)}$ whose vertex set is the set of of all 
$(k+1)$- cliques in $G$ and in which two vertices (i.e., $(k+1)$-cliques) are \emph{adjacent} when they share a common face that has a minimum of $l$ vertices, where $k \geq 1$ and $1 \leq l \leq k$. Such an adjacency matrix is built for each $k$-skeleton and therefore $\mathcal{X}(G)$ is expressed as a set of matrices $\{G^{(k,l)}  \}_{k=0}^h$, where $(k+1)$ is the dimension of the cliques.


\begin{algorithm}[tbp]
	\caption{Matching Random Clique Complexes}
	\label{alg:matching}
	\begin{flushleft}
		\textbf{Input:} $\mathcal{X}(G) = \{  G^{(k,l)}   \}_{k=0}^h$ and 
		$\mathcal{X}(G') = \{ G'^{(k,l)}   \}_{k=0}^h$
	\end{flushleft}
	\begin{algorithmic}[1]
	
		\FOR{ $k = h \dots 0$}
		\STATE Let $M, M'$ be the total number of $(k+1)$-cliques in $G^{(k,l)}$ 
		and $G'^{(k,l)}$, respectively
	
		\STATE $\mathcal{L} := \{   c_i^{(k)}  \}_{i=0}^{M-1} $  \COMMENT {list of barycenters}
	
		\FOR{ $i = 0 \dots M-1$} 
		\STATE $\mathcal{N}_i := \mathcal{N}_i \cup \left\{  g^{(k,l)}_{(x,:)} \mid x=i,   g^{(k,l)}_{(x,y)} \neq 0   \right\}  $
		\STATE $\mathcal{N} := \mathcal{N} \cup \{  \mathcal{N}_i \} $ \COMMENT {clique neighborhoods}
		\ENDFOR
	
		\FOR{ $i = 0 \dots M-1$} 
		\STATE $\alpha_i := [ \alpha_1, \dots, \alpha_{M-1}]^T$
		\STATE $\alpha := \alpha \cup \{ \alpha_i \}$  \COMMENT{affine weight vectors}
		\ENDFOR
	
		\STATE Repeat steps $3$--$11$ on $G'^{(k,l)}$ for $\mathcal{L}',\mathcal{N}'$ and $\alpha'$.
	
		\STATE Build cost matrix $C^{(k)}$ from weights vectors $\alpha, \alpha'$
		\STATE $X^{*}_k := $ \textbf{Kuhn-Munkres }($G^{(k,l)},G'^{(k,l)},C^{(k)} $)	
		\ENDFOR
	\end{algorithmic}
	\begin{flushleft}
		\textbf{Return:}  $ \{ X^{*}_0, \dots, X^{*}_h \}$ \text{ \# set of permutation matrices}
	\end{flushleft}
\end{algorithm}
	
\subsection{Problem Setup}
\label{subsec:setup}
The problem of matching random clique complexes each of dimension $h$, is the estimation of a set of 
optimal bijective maps of the form $\mathcal{M}_i: \mathcal{X}^{(i)}(G) \rightarrow \mathcal{X}^{(i)}(G')$, for all $i \leq h$, subject to assignment constraints.
This can be formulated as a \emph{constrained quadratic assignment problem}, which can later be relaxed to a linear programming optimization problem.

Given two $h$-dimensional random clique complexes $\mathcal{X}(G) = \{  G^{(k,l)}   \}_{k=0}^h$ and 
$\mathcal{X}(G') = \{ G'^{(k,l)}   \}_{k=0}^h$, let $X=\{  X_0, \dots, X_h   \} \in \Pi$ be a set of permutation matrices
such that $X_k$ encodes assignments/matchings from $G^{(k,l)}$ to $G'^{(k,l)}$. 
The combinatorial matching requires the optimal set of permutation matrices that best \emph{align} 
$\mathcal{X}(G)$ and $\mathcal{X}(G')$. More formally, this can be expressed as the following constrained optimization problem
\begin{equation}
\label{eq:opt}
\begin{aligned}
& \operatornamewithlimits{argmin}_{X_0,\dots,X_h}
& &   \sum_{k=0}^{h} \lVert  G^{(k,l)} X_k - X_k G'^{(k,l)}    \rVert_F^2 \\
& \text{subject to}
& & \forall k \leq h, \mathds{1}^T X_k = \mathds{1},  X_k^T \mathds{1} = \mathds{1}
\end{aligned}
\end{equation}

\subsection{Our Algorithm}
\label{subsec:algo}

At a high level, our goal is to minimize $\lVert \mathcal{X}(G) - \mathcal{X}(G') \rVert_{\mathcal{C}}$, where $\mathcal{C}$ is a combinatorial distance between two random clique complexes. Traditional metrics like Hausdorff distance are not suitable here because random clique complexes are combinatorial topological spaces.
Recall that $\mathcal{X}(G)$ is comprised of a family of $k$-skeletons $\{  \mathcal{X}^{(k)}(G)  \}_{k=0}^h$, where each $k$-skeleton contains cliques whose dimension is at most $k+1$ and $\mathcal{X}^{(k)}(G)$ has a maximum dimension $h$.
The solution of the optimization problem outlined in Equation~(\ref{eq:opt}) aims to find a set of permutation matrices $\{X_1, \dots, X_h \}$ that minimizes the overall number of \emph{misalignments} between equi-dimensional faces of $\mathcal{X}(G)$ and $\mathcal{X}(G')$, i.e., cliques belonging to the corresponding $k$-skeletons, and thus producing the optimal least cost assignment between $\mathcal{X}(G)$ and $\mathcal{X}(G')$. 

Algorithm~\ref{alg:matching} presents our method to solve the combinatorial optimization problem (Equation~(\ref{eq:opt})). In decreasing order of clique dimensionality, for a fixed dimension $k$ and given the adjacency matrices $G^{(k,l)}$ and $G'^{(k,l)}$ for $k$-skeletons $\mathcal{X}^{(k)}(G)$ and $\mathcal{X}^{(k)}(G')$, respectively. In every iteration, our objective is to solve $\operatornamewithlimits{argmin}_{X_k} \lVert G^{(k,l)}X_k - X_k G'^{(k,l)}  \rVert_F^2$ to find the optimal permutation $X^{*}_k$.
We assume the barycenters of every clique is pre-computed (Step~$3$). 
Next, the \emph{neighborhood} $\mathcal{N}_i$ of the $i$-th clique is computed as the set of entries with $1$s in the $i$-th row of $G^{(k,l)}$  (Step~$5$). We denote the collection of every clique's neighborhood as $\mathcal{N}$ (Step~$6$). 
An important objective of our method is to capture the \emph{geometric properties} of the neighborhood of every clique. We achieve this by characterizing the $i$-th clique's barycenter $c_i^{(k)}$ as an \emph{affine combination} of the barycenters (in all dimensions) associated with the cliques in its corresponding neighborhood $\mathcal{N}_i$.
Given an arbitrary clique's barycenter $c_i^{(k)}$, let $\{ x_1^{(k)}  , \dots,  x_n^{(k)} \}$ denote the barycenters of its $n$ adjacent cliques. Then, $c_i^{(k)}$ expressed as $\sum_{i=1}^{n} \alpha_i x_i^{(k)}$ is an affine combination of the $x_i^{(k)}$s, if $\sum_{i=0}^{n} \alpha_i = 1$, i.e., the weights $\alpha_i$ sum to $1$.
Among all possible affine representations of $c_i^{(k)}$ we chose to use \emph{least squares} to guarantee minimal error under L$2$-norm, and furthermore it assigns non-zero weights to each of its adjacent clique barycenters, thereby capturing the local geometric properties in its neighborhood. 
The weight vector $\alpha_i$ is then calculated for each clique  (Step~$9$) and $\alpha$ denotes a collection of such weight vectors (Step~$10$). Next, a cost matrix is built by computing the L$2$-norm distance between weight vectors $\alpha$ and $\alpha'$ (Step~$13$). Finally, the \emph{Kuhn-Munkres}~\cite{Kuhn1955} algorithm is invoked with both the adjacency matrices and the cost matrix, which arrives at the optimal assignment (Step~$14$ ). At the end of all iterations, our method returns a set of optimal assignments for matches between each $k$-skeleton for every dimension below $h$ and the algorithm terminates. We refer the reader to our supplementary section for a working example.

\subsection{Complexity Analysis}
\label{subsec:complexity}

To begin our analysis, we must first ascertain the dimensionality of $G^{(k,l)}$, which is governed by the total number of $(k+1)$-cliques that exist in the underlying random graph $G$. It is important to note that there doesn't exist any closed form solution to counting the number of cliques of a given dimension in $G$.

We consider the distribution of a random variable $X_n(k)$ counting the number of $(k+1)$-cliques in a realization of $G$. We show in Appendix A of our supplementary material that this count is upper bounded by 
$(en/k)^k$, where $e$ is \emph{Euler's number}. 
This can be expressed in asymptotic notation as $O(n^k)$.
As dimensionality increases, there occurs an explosion in the number of cliques.
Fortunately, $G^{(k,l)}$ is a sparse matrix and its \emph{effective} dimensionality measured by the number of non-zero rows, i.e., the number of cliques with non-empty neighborhoods, is of order $O(nnz(G^{(k,l)}))$. Therefore, we set out to count the number of non-zero entries in $G^{(k,l)}$.

We use a seminal result by Bollob{\'{a}}s~\cite{Bollobas2009}, where they identify a threshold probability for \emph{percolation} of cliques in $G$ for all fixed $k$ and $l$, which is given by 
$p = \Theta \left(n ^{\frac{-2}{k+l-1}} \right)$. Moreover, they proved that for $p$ around this threshold, the number of cliques asymptotically converge to a \emph{Poisson} distribution. Exceeding this threshold results in formation of \emph{giant connected clique clusters}, which  causes an explosion in the number of possible cliques. 

Recall from our definition of $G^{(k,l)}$, that two cliques are adjacent if they share at least $l$ vertices. In order to analyze this further, we imagine an entry in $G^{(k,l)}$ occurs when we can migrate a $(k+1)$-clique from its original position to an adjacent clique by relocating exactly $(k+1-l)$ vertices and leaving the remaining $l$ vertices intact. The expected number of such relocations is given by $\left( {k+1 \choose l} -1 \right)  {n \choose k+1-l}  p^{  \left(  {k+1 \choose 2} - {l \choose 2}  \right)}$,
where the first term denotes the number of possible vertices in a $(k+1)$-clique 
that can be chosen for relocation, the second term counts the number of new adjacent positions a clique can relocate to, and the final term decides the probability of relocations that are correct and acceptable.
In our case, we define cliques to be adjacent to one another when they share at least $l=k$ vertices. This is done in order to keep the number of adjacent cliques to a manageable size during experiments.
Setting $l=k$, gives $kn p^k$ expected relocations, which in turn estimates $nnz(G^{(k,l)})$.


Note that for every iteration in Algorithm~\ref{alg:matching}, the dominating cost is that of running the Kuhn-Munkres matching algorithm in Step~$14$, which has a \emph{cubic} cost in $nnz(G^{(k,l)})$.
Let $\mathcal{C}_{max}$ denote an upper bound on all the number of non-zero entries in $ \{  G^{(k,l)}  \}_{i=0}^h $.
Then, every iteration has a runtime $O(\mathcal{C}_{max}^3)$ and therefore after $h$ iterations the final cost is $O(h\mathcal{C}_{max}^3)$. Observe that as the dimensionality of the cliques increases in every iteration, 
$p^k$ decays very sharply and hence drastically reduces $nnz(G^{(k,l)})$, which in turn reduces the overall matching cost. Finally, the storage complexity can simply be given as $O(h\mathcal{C}_{max})$.

\section{Theoretical Analysis of QAP}
\label{sec:theory}
In this section, we present three related results in the context of matching random matrices, namely: (i) concentration inequality of eigenvalue bounds on the \emph{QAP trace formulation} for random symmetric matrices, (ii) tighter concentration inequality of eigenvalue bounds on the \emph{Lawler QAP formulation} on random symmetric matrices in the context of works that use \emph{affinity matrices}, and (iii) provide an asymptotic analysis on the worst to best case ratio of a QAP for higher-dimensional clique adjacency matrices. For ease of notation, we will refer to the random clique adjacency matrices simply as $A$ and $B$.

\subsection{Eigenvalue Bounds of Trace QAP Formulation on Random Matrices}

Let $A=(a_{vv'})$, $B=(b_{vv'}) \in \mathds{R}^{n \times n}$ be random real-symmetric matrices. Let $X=(x_{ij}) \in \mathds{R}^{n \times n}$ be a \emph{permutation matrix}. 
Then, the trace formulation of a QAP is given by
\[
\text{minimize } \tr(AXBX^T) \text{ s.t. } X \in \Pi_X
\]
where $\Pi_X$ is the set of permutation matrices.

Let 
$\lambda_1(A) \leq \lambda_2(A) \leq \dots \leq \lambda_n(A)$ and 
$\lambda_1(B) \geq \lambda_2(B) \geq \dots \geq \lambda_n(B)$\footnote{The two sets of eigenvalues differ in ordering.}
be the eigenvalues of $A$ and $B$, respectively. 
Let the corresponding eigen-decompositions of matrices $A$ and $B$, be given by 
$A=Q_A \Lambda_A Q_A^T$ and $B=Q_B \Lambda_B Q_B^T$, where 
$\Lambda_A = \diag(\lambda_1(A),\dots,\lambda_n(A))$ and 
$\Lambda_B = \diag(\lambda_1(B),\dots,\lambda_n(B))$ with their corresponding 
orthogonal eigenvector matrices $Q_A$ and $Q_B$.
Finke et. al.~\cite{martello1987} gave the following eigenvalue bounds.
	\begin{theorem}
	\label{thm:1}
	Let $A$ and $B$ be symmetric matrices. Then for all $X \in X_\Pi$, 
	\begin{enumerate}
		\item 	$\tr(AXBX^T) = \lambda(A)^T Q^{(X)} \lambda(B)$, \\
		where $Q^{(X)} = \langle Q_A^{(i)} , X Q_B^{(j)} \rangle^2 $ with 
		vectors of eigenvalues given by $\lambda(A) = (\lambda_i(A))$ and $\lambda(B) = (\lambda_i(B))$. $Q_A^{(i)}$ and $Q_B^{(i)}$ denote the $i$-th eigenvectors of
		$A$ and $B$, respectively;
		\item 
		$ \mathcal{L} \leq \tr(AXBX^T) \leq \mathcal{U}$, where \\
		$\mathcal{L}=\sum_{i=1}^{n} \lambda_i(A) \lambda_i(B)$, and\\
		$\mathcal{U}=\sum_{i=1}^{n} \lambda_{n-i+1}(A) \lambda_i(B)$
	
	
	
	\end{enumerate}
\end{theorem}

\iffalse
Given $\tr(AXBX^T) = \lambda(A)^T Q^{(X)} \lambda(B)$, 
where $Q^{(X)} = \langle Q_A^{(i)} , X Q_B^{(j)} \rangle^2 $ with 
vectors of eigenvalues given by $\lambda(A) = (\lambda_i(A))$ and $\lambda(B) = (\lambda_i(B))$. Furthermore, $Q_A^{(i)}$ $Q_B^{(i)}$ denote the $i$-th eigenvectors of
$A$ and $B$, respectively; then 
\begin{align*}
\sum_{i=1}^{n} \lambda_i(A) \lambda_i(B) 
\leq \tr(AXBX^T)
\leq \sum_{i=1}^{n} \lambda_{n-i+1}(A) \lambda_i(B) 
\end{align*}
\fi

It was also noticed by Finke~\cite{martello1987} that these bounds can further be 
tightened by reducing the \emph{spreads} of matrices $A$ and $B$, where the spread of a matrix $A$, denoted by $\mathfrak{S}(A)$, is given by $\mathfrak{S}(A) = \max_i {\lambda_i(A)} - \min_i {\lambda_i(A)}$.
There is no formula to compute the spread of a matrix directly, 
so Finke et. al.~\cite{martello1987} suggested a 
\emph{reduction} method to further sharpen the bound by replacing matrices $A$ and $B$ by \emph{smaller spread} symmetric matrices $\tilde{A}$ and $\tilde{B}$.
The reductions are achieved as
$\tilde{A} = A - M_A - M_A^T - \mathfrak{D}_A$ and 
$\tilde{B} = B - M_B - M_B^T - \mathfrak{D}_B$  	
, where $M_A$, $M_B$ are matrices with constant columns and 
$\mathfrak{D}_A$, $\mathfrak{D}_B$ are diagonal matrices, whose values 
are chosen appropriately in order to tighten the bounds on spreads $\mathfrak{S}(\tilde{A})$
and $\mathfrak{S}(\tilde{B})$.

\textbf{Our bounds on Random Matrices}: We propose new measure concentration 
inequalities on the spread of a random matrix, by redefining the spread in an alternate fashion that is more amenable to our analysis. 
Consider our reduced random symmetric matrix $\tilde{A} \in \mathds{R}^{n \times n}$, with eigenvalues 
$\lambda_1(\tilde{A}) \leq \dots \leq \lambda_n(\tilde{A})$, we define the \emph{gap} (spacing) between its consecutive eigenvalues as $\delta_i(\tilde{A}) := |\lambda_{i+1}(\tilde{A}) - \lambda_i(\tilde{A})|$ for $1 \leq i \leq n-1$. Then, the spread $\mathfrak{S}(\tilde{A})$ for the reduced matrix $\tilde{A}$ can be redefined as:
$\mathfrak{S}(\tilde{A}) = \sum_{i=1}^{n} \delta_i(\tilde{A})$

We begin by upper bounding $\delta_i(\tilde{A})$ using the following lemma~\ref{lemma:2} (proof in supplementary notes).
\begin{lemma}
	\label{lemma:2}
	Let $\vertiii{.}$ denote an algebraic matrix norm on a space of real $n \times n$ matrices $\mathcal{M}_n$, then for any $A \in \mathcal{M}_n$, 
	\begin{align*}
	\delta_i(A) \leq 2 \vertiii{A} 
	\end{align*}
\end{lemma}
\iffalse
	\begin{proof}
	For ease of understanding, we drop the $(A)$ as it is obvious from context. Let $\lambda_i$ be the $i$-th eigenvalue of $A$, and let $x_i \neq 0$ be its corresponding eigenvector. From $A x_i = \lambda_i x_i$, we have
	\begin{align*}
	A X_i = \lambda_i X_i, \text{\hspace{1em}where } 
	X_i:= \begin{bmatrix}[c|c|c]
	x_i & \dots & x_i
	\end{bmatrix}
	\in \mathcal{M}_n \setminus \{0\}
	\end{align*} 
	It follows,
	\begin{align*}
	|\lambda_i|  \vertiii{X_i} \vertiii{X_i}  = \vertiii{\lambda_i X_i}    
	=  \vertiii{AX_i} \leq \vertiii{A} \vertiii{X_i} .
	\end{align*} 
	
	As $\vertiii{X}$ is non-negative, we get $|\lambda_i| \leq \vertiii{A}$.
	Thus, every eigenvalue of $A$ is upper bounded by the matrix norm $\vertiii{A}$.
	Applying the triangle inequality, we get that $\delta_i(A) = |\lambda_{i+1}(A) - \lambda_i(A)| \leq 2 \vertiii{A}$, which completes the proof.
\end{proof}

\begin{proof}
	For ease of understanding, we drop the $(A)$ as it is obvious from context. Let $\lambda_i$ be the $i$-th eigenvalue of $A$, and let $x_i \neq 0$ be its corresponding eigenvector. From $A x_i = \lambda_i x_i$, we have
	\begin{align*}
	A X_i = \lambda_i X_i, \text{\hspace{3em}where } 
	X_i:= \begin{bmatrix}[c|c|c]
	x_i & \dots & x_i
	\end{bmatrix}
	\in \mathcal{M}_n \setminus \{0\}
	\end{align*} 
	It follows,
	\begin{align*}
	|\lambda_i|  \vertiii{X_i} \vertiii{X_i}  = \vertiii{\lambda_i X_i}    
	=  \vertiii{AX_i} \leq \vertiii{A} \vertiii{X_i} .
	\end{align*} 
	
	As $\vertiii{X}$ is non-negative, we get $|\lambda_i| \leq \vertiii{A}$.
	Thus, every eigenvalue of $A$ is upper bounded by the matrix norm $\vertiii{A}$.
	Applying the triangle inequality, we get that $\delta_i(A) = |\lambda_{i+1}(A) - \lambda_i(A)| \leq 2 \vertiii{A}$, which completes the proof.
\end{proof}
\fi

To the best of the author's knowledge there does not exist a known distribution of eigenvalue gaps for a symmetric random matrix. We now attempt to give concentration inequalities for the tail probabilities of the sum of eigenvalue gaps, i.e., the spread.
For our i.i.d. random matrix $A \in \mathds{R}^{n \times n}$, 
consider the sequence of independent eigenvalue gaps 
$\delta_1(A), \dots ,\delta_n(A)$, where each $\delta_i(A)$ is upper bounded by $2 \vertiii{A}$,
as shown in Lemma~\ref{lemma:2}. Let us denote their sum as 
$\mathcal{S}_n(A) := \delta_1(A) + \dots + \delta_n(A)$.
As $\delta_1(A), \dots ,\delta_n(A)$ are independent scalar random variables with 
$\delta_i(A) \leq \vertiii{A}$ a.s, with mean $\mu_i(A)$ and variance
$\sigma_i^2(A)$. Then, using Chernoff bounds, for any $ \epsilon>0$, we have
\begin{align*}
\mathbb{P}(| \mathcal{S}_n(A) - \mu | \geq \epsilon \sigma) 
\leq K  \max \left( e^{(-p \epsilon^2  ) } , e^{  ( -p \epsilon \sigma / 2 \vertiii{A} ) }  \right)
\end{align*} 
for some absolute constants $K,p >0$.
The Chernoff inequality above, shows that $\mathcal{S}_n(A)$ is sharply concentrated in the range
$n \mu + O(\sigma\sqrt{n})$, when $\epsilon$ is not too large.

\subsection{Eigenvalue bounds on Lawler's QAP on Random Affinity Matrices}
In literature, many graph matching algorithms use Lawler's QAP formulation. 
Recall, $A=(a_{ij}) ,B=(b_{uv}) \in \mathds{R}^{n \times n}$. 
Let $\Omega(a_{i,j}, b_{u,v})$ denote the \emph{pairwise affinity score} of 
assigning the $(i,j)$-th entry in $A$ to the $(u,v)$-th entry in $B$, implying that 
node $a_i$ is matched to node $b_u$ and node $a_j$ to node $b_v$, simultaneously.
Then, the \emph{affinity matrix} $\mathscr{A} \in \mathds{R}^{n^2 \times n^2}$ is given by 
$\mathscr{A}[(i-1)n + u, ( j-1  )n + v]=  \Omega(a_{i,j}, b_{u,v})$ and the optimal assignment to Lawler's QAP is the one that maximizes the sum total pairwise affinity scores. 
Leordeanu et. al.~\cite{Leordeanu2005} show via a spectral relaxation 
that Lawler's WAP reduces to solving 
$w^* = \argmaxI_w \frac{w^T Aw}{w^Tw} $, $w \in \mathds{R}^{n^2}$. This is solved by finding the leading eigenvalue $\lambda_1(\mathscr{A} )$.

As illustrated in~\cite{Alon2002}, we also use \emph{Talagrand's concentration inequality}~\cite{Talagrand1995}. We provide a tighter bound in the case of our affinity matrix using \emph{Rayleigh's quotient}. 
\begin{theorem}
	\label{thm:bound}
	For a random affinity matrix $\mathscr{A} \in R^{m \times m}$ and for a positive constant $t$,

	$\mathbb{P}[ |\lambda_1( \mathscr{A}) - \mathcal{M}| \geq t] \leq 4e^{-t^2/8}$,

	where $\mathcal{M}$ is the median of $\lambda_1( \mathscr{A})$.\qed
\end{theorem}

\textbf{Discussion:} We further investigate the robustness of affinity-matrix based graph matching solutions when dealing with \emph{missing} or \emph{incomplete} data. We show the sharpness of our result in Theorem~\ref{thm:bound} on the affinity matrix, similar to~\cite{Alon2002}, who analyze their results using \emph{fat} matrices as an example. 
Consider an affinity matrix $\mathscr{A}=(a_{ij}) \in \mathcal{R}^{m \times m}$, whose entries are i.i.d. 
Bernoulli distributed. Simulating missing affinity scores due to missing edge assignments, we set $a_{ij}=1$ with probability $1/4$ and $a_{ij}=0$ with probability $3/4$. Notice that our random affinity matrix represents the Erdős-Rényi graph $G(m,1/4)$.  
As shown in~\cite{Alon2002}, the median and expectation of $\lambda_1(A)$ differ by a constant factor. Let $G=(V,E)$ denote a general undirected graph, where the degree of each vertex $v \in V$ is given by $d_G: V \rightarrow \mathds{Z}$. Then, the \emph{average degree} of $G$ is given by $\bar{d} = \sum_{v \in V} d_G(v) / |V|$ and its \emph{maximum degree} is 
$\Delta = \max_{v \in V} d_G(v)$. It is then well known that $\bar{d} \leq \lambda_1(A) \leq \Delta$, i.e., the largest eigenvalue of a graph is squeezed between its average and maximum degree. 

Let $|E|$ denote the total number of edges in $G(m,1/4)$, then the \emph{average degree} of $G(m,1/4)$ is given by $2 |E|/n$, where $ |E| = (Bin( {m \choose 2}, 1/4 )$ and the standard deviation of $|E|$ is $ \sqrt{ {m \choose 2} (1/4) (3/4) } = \Theta(m)$. 
For large $m$, our binomial distribution converges to a normal distribution. Therefore, we calculate the probability for the total number of edges $|E|$ to deviate from its expectation by $t$ standard deviations as $e^{-\Theta(t^2)}$. Furthermore, we know that if $|E|$ exceeds its expectation by $\Theta(tn)$, then the average degree $\bar{d}$ must also correspondingly exceed its expectation by $\Theta(t)$. Therefore, the probability of the average degree $\bar{d}$ exceeding its expectation by $t$ standard deviations is \emph{at least} $e^{- \Theta(t^2)}$. Given that $\bar{d} \leq \lambda_1(A)$, it follows that $\lambda_1$ exceeding its expectation is also lower bounded by the same $e^{- \Theta(t^2)}$. The bounds achieved are \emph{tight up to a constant factor} in the exponent. Our experimental results on \emph{Factorized Graph Matching (FGM)} by Zhou et. al.~\cite{ZhouD16} and \emph{Re-weighted Random Walk Matching (RRWM)}~\cite{cho2010reweighted}  also support the finding that affinity matrix based matching solutions are more robust to missing edges due to occlusions in data. 

\subsection{Asymptotic Analysis of Higher-Order Clique Assignment}	
Following along the same lines as Finke et. al.~\cite{martello1987}, we study the asymptotic behavior of the worst to most optimal ratio and present it as the following theorem. 
\begin{theorem}
	Given random clique adjacency matrices $A^{k,l}_n, B^{k,l}_n \sim Pois(\lambda) $ and their associated cost matrix $\mathcal{C}_{vv'} \sim Pois(\lambda)$. We denote by $\lambda_e := \mathbb{E}(\mathcal{C}_{vv'}) $ and 
	$\lambda_v := \mathbb{V}(\mathcal{C}_{vv'})$ the expectation and variance of our Poisson distributed cost function. 
	For $\epsilon>0$ and $p = \Theta \left( n^{\frac{-2}{(k+l-1)}} \right)$, we have the following bound on the ratio of the worse to the best solution as
	\begin{align*}
	&\mathbb{P}\left\lbrace 
	\frac{\max\limits_{\pi \in \Pi}  \sum\limits_{vv'} \mathcal{C}_{vv'}}  
	{\min\limits_{\pi \in \Pi}  \sum\limits_{vv'} \mathcal{C}_{vv'}}
	\leq 1 + \epsilon	 
	\right\rbrace \\
	&\geq 1 - 
	2|\Pi|\exp \left( -2|S_\pi| \left( \frac{\epsilon'\sqrt{\lambda_v}}
	{\epsilon' + 2\lambda_v}   \right)^2 \right) =: \psi(n,\epsilon)
	\end{align*}
	where, $|\Pi| = n!$, $|S_\pi|= {{n+1} \choose 2}$, $\lim_{n \to \infty} \psi(n,\epsilon) = 1$\qed
\end{theorem}

\begin{figure*}[t!]
	\centering
	\subfigure[]{%
		\label{fig:zero11}%
		\includegraphics[width=27mm]{Experiments/Transformation/house0.png}}%
	\qquad
	\subfigure[]{%
		\label{fig:first11}%
		\includegraphics[width=27mm]{Experiments/Transformation/HouseRotate180.png}}%
	\qquad
	\subfigure[]{%
		\label{fig:second11}%
		\includegraphics[width=27mm]{Experiments/Transformation/HouseMirrorx_marked.png}}%
	\qquad
	\subfigure[]{%
		\label{fig:third11}%
		\includegraphics[width=27mm,height=27mm]{Experiments/Transformation/HouseScale05_15.png}}%
	\qquad
	\subfigure[]{%
		\label{fig:fourth11}%
		\includegraphics[width=27mm,height=27mm]{Experiments/Transformation/HouseShear05_marked.png}}%
	\vspace{-5mm}
	\caption{(a) Original House Frame, (b)--(e) Four transformations on house frame: Rotation, Reflection, Scaling and Shear (green markers show true matching case with the original frame (a)).}
	\label{fig:HouseTransform}
\end{figure*}
\begin{table*}[h]
	{\tiny
		\hfill{}
		\begin{tabular}{l||l|l|l|l|l|l|l|l|l|l}
			\hline
			\multicolumn{1}{l||}{\textbf{Algorithms}} & \multicolumn{2}{|l|}{\textbf{$20^{\circ}$} Rotation} & \multicolumn{2}{|l|}{\textbf{$60^{\circ}$} Rotation}&  \multicolumn{2}{|l|}{\textbf{Reflection}} & \multicolumn{2}{|l|}{\textbf{Scaling}} & \multicolumn{2}{|l}{\textbf{Shear}}\\
			
		
			\cline{1-11}
			& 20\% & 40\% & 20\% & 40\% & 20\% & 40\% & 20\% & 40\% & 20\% & 40\% \\
			\cline{2-11}
			\hline
			\hline
			OurMethod & \textbf{0.01 $\pm$ 0.0} & \textbf{0.01 $\pm$ 0.0} & \textbf{0.05 $\pm$ 0.0} & \textbf{0.03 $\pm$ 0.0} & \textbf{0.0 $\pm$ 0.0} & \textbf{0.0 $\pm$ 0.0} & \textbf{0.1 $\pm$ 0.0} & \textbf{0.2 $\pm$ 0.1} & \textbf{0.0 $\pm$ 0.0} & \textbf{0.0 $\pm$ 0.0}\\
			EigenAlign & 63.97 $\pm$ 1.0 & 70.19 $\pm$ 1.0 & 65.39 $\pm$ 1.5 & 70.88 $\pm$ 1.9 & 62.5 $\pm$ 0.4 & 64.97 $\pm$ 0.2 & 62.37 $\pm$ 0.6 & 65.42 $\pm$ 1.1 & 61.04 $\pm$ 0.2 & 62.36 $\pm$ 0.7\\
		
			FGM & 3.6 $\pm$ 0.5 & 7.4 $\pm$ 0.5 & 18.0 $\pm$ 0.0 & 36.4 $\pm$ 0.5 & \textbf{0.0 $\pm$ 0.0} & \textbf{0.0 $\pm$ 0.0} & 2.20 $\pm$ 1.0 & 3.4 $\pm$ 0.5 & \textbf{0.0 $\pm$ 0.0} & \textbf{0.0 $\pm$ 0.0}\\
			LAI-LP & 40.90 $\pm$ 0.9 & 43.07 $\pm$ 1.2 & 49.9 $\pm$ 1.2 & 61.16 $\pm$ 0.7 & 49.07 $\pm$ 0.9 & 59.41 $\pm$ 0.8 & 43.12 $\pm$ 5.5 & 45.91 $\pm$ 3.2 & 39.04 $\pm$ 0.7 & 38.27 $\pm$ 0.6\\
			PermSync & 13.36 $\pm$ 1.1 & 16.48 $\pm$ 0.4 & 25.48 $\pm$ 1.4 & 41.59 $\pm$ 0.7 & 26.24 $\pm$ 3.4 & 42.31 $\pm$ 5.1 & 12.15 $\pm$ 0.7 & 13.7 $\pm$ 0.6 & 10.19 $\pm$ 0.8 & 6.23 $\pm$ 3.2\\
		
			RRWM & 2.0 $\pm$ 0.0 & 4.0 $\pm$ 0.0 & 14.0 $\pm$ 0.0 & 27.0 $\pm$ 0.0 & \textbf{0.0 $\pm$ 0.0} & \textbf{0.0 $\pm$ 0.0} & 5.0 $\pm$ 2.0 & 10.0 $\pm$ 1.5 & \textbf{0.0 $\pm$ 0.0} & \textbf{0.0 $\pm$ 0.0}\\
			Tensor & 6.88 $\pm$ 0.0 & 14.39 $\pm$ 0.8 & 19.24 $\pm$ 0.6 & 37.95 $\pm$ 0.8 & 17.1 $\pm$ 0.5 & 34.69 $\pm$ 0.4 & 2.91 $\pm$ 0.6 & 6.19 $\pm$ 0.7 & \textbf{0.0 $\pm$ 0.0} & \textbf{0.0 $\pm$ 0.0}\\
			IPFP & 6.8 $\pm$ 1.5 & 11.6 $\pm$ 0.5 & 18.2 $\pm$ 0.5 & 35.0 $\pm$ 0.0 & 1.0 $\pm$ 0.0 & 1.0 $\pm$ 0.0 & 7.4 $\pm$ 2.0 & 15.6 $\pm$ 1.5 & 0.8 $\pm$ 0.5 & 0.6 $\pm$ 0.5\\
			PM & 43.2 $\pm$ 1.0 & 49.4 $\pm$ 1.5 & 46.4 $\pm$ 1.0 & 56.8 $\pm$ 0.5 & 37.0 $\pm$ 0.0 & 37.0 $\pm$ 0.0 & 45.6 $\pm$ 0.5 & 53.4 $\pm$ 0.5 & 37.8 $\pm$ 0.5 & 38.4 $\pm$ 0.5\\
		
			SMAC & 12.6 $\pm$ 0.5 & 20.0 $\pm$ 0.0 & 19.0 $\pm$ 0.0 & 33.0 $\pm$ 0.0 & 4.0 $\pm$ 0.0 & 4.0 $\pm$ 0.0 & 13.8 $\pm$ 2.5 & 23.4 $\pm$ 0.5 & 4.0 $\pm$ 0.0 & 4.2 $\pm$ 0.5\\
			SM & 32.6 $\pm$ 0.5 & 39.4 $\pm$ 1.5 & 37.8 $\pm$ 1.0 & 49.6 $\pm$ 1.0 & 25.0 $\pm$ 0.0 & 25.0 $\pm$ 0.0 & 33.0 $\pm$ 1.5 & 40.4 $\pm$ 0.5 & 25.2 $\pm$ 0.5 & 24.8 $\pm$ 0.5\\
			GA & 34.0 $\pm$ 1.5 & 37.0 $\pm$ 1.0 & 38.4 $\pm$ 1.0 & 47.0 $\pm$ 0.0 & 30.0 $\pm$ 0.0 & 30.0 $\pm$ 0.0 & 37.6 $\pm$ 1.0 & 45.4 $\pm$ 0.5 & 29.0 $\pm$ 0.0 & 28.4 $\pm$ 0.5\\
			Munkres & 35.25 $\pm$ 1.5 & 36.47 $\pm$ 0.9 & 44.75 $\pm$ 0.8 & 55.13 $\pm$ 1.3 & 46.05 $\pm$ 1.2 & 57.25 $\pm$ 2.2 & 37.32 $\pm$ 1.5 & 37.81 $\pm$ 1.8 & 32.84 $\pm$ 0.3 & 31.27 $\pm$ 0.8\\			
			\hline
	\end{tabular}}
	\hfill{}
	\caption{Error (\%) of transformation on CMU House: inserted $20\%$ and $40\%$ impurity in CMU House frame sequence randomly for rotation ($20^{\circ}$, $60^{\circ}$), reflection, scaling and shear. Minimum error (\%) is shown in bold. Matching is computed for $111$ frames from the $1^{st}$ frame to the other $110$ frames. Our method shows best performance among all the methods.}
	\label{tb:Transformation}
\end{table*}


\section{Experiments}
\label{sec:exp}
Here, we study the robustness of various matching algorithms when affected by missing or incomplete information and transformations (both affine and non-affine) on synthetic and real-world datasets. For the sake of brevity, we report detailed dataset descriptions in our supplementary notes.
The graph matching algorithms can broadly be classified based on their use of
(i) \emph{affinity-matrix:} FGM~\cite{ZhouD16,zhou2013deformable}\footnote{\href{http://www.f-zhou.com/gm\_code.html}{FGM}}, RRWM~\cite{cho2010reweighted}, 
(ii) \emph{Eigenvalues:} EigenAlign~\cite{feizi2016spectral}\footnote{\href{https://github.com/SoheilFeizi/spectral-graph-alignment}{EigenAlign}}, SM~\cite{Leordeanu2005}, SMAC~\cite{cour2007balanced}, PermSync~\cite{pachauri2013solving}\footnote{\href{http://pages.cs.wisc.edu/~pachauri/perm-sync/}{PermSync}},
(iii) \emph{LP relaxation:} GA~\cite{gold1996graduated}, Kuhn-Munkres~\cite{Kuhn1955}, 
(iv) \emph{Integer QAP:} IPFP~\cite{leordeanu2009integer}, 
(v) \emph{Probabilistic matching:} PM~\cite{zass2008probabilistic}, 
(vi) \emph{Higher-order matching given complete data:} Tensor~\cite{duchenne2011tensor}\footnote{\href{http://www.cs.cmu.edu/~olivierd/}{Tensor}},
and (vii) \emph{Geometric and Feature matching:} LAI-LP~\cite{li2013object}\footnote{This algorithm serves as our \emph{naive baseline} as it directly uses neighborhood properties of the underlying graph (\href{http://www.ee.cuhk.edu.hk/~hsli/}{LAI-LP}).}. Our code\footnote{\href{https://github.com/charusharma1991/RandomCliqueComplexes_ICML2018}{Our Method}} is publicly available.


\subsection{Effect of Affine Transformations}
\textbf{Simulated Dataset:} We perform affine transformations on \href{http://vasc.ri.cmu.edu/idb/html/motion/house/index.html}{CMU House}, which is a sequence of $N$ frames extracted from a video. More specifically, we uniformly sample frames (at $20\%$ and $40\%$) and perform affine transformations on the selected frames to distort them. Figure~\ref{fig:HouseTransform} shows examples of affine transformations on house frame sequences. Table \ref{tb:Transformation} shows the comparative error in matching for all the algorithms. We now describe each affine transformation as performance metrics in our experiments.
\begin{center}
	\begin{table*}
		{\tiny
			\hfill{}
			\begin{tabular}{l||l|l|l|l|l|l}
				\hline
				\textbf{Methods} & \textbf{Car} & \textbf{Bike} & \textbf{Butterfly} & \textbf{Magazine} & \textbf{Building} & \textbf{Book}\\
				\hline
				\hline
				OurMethod & \textbf{4.14 $\pm$ 2.45}/ 7.13 & \textbf{3.15 $\pm$ 0.32}/ 6.96 & \textbf{3.89 $\pm$ 0.23}/ 14.76 & \textbf{0.48 $\pm$ 0.02}/ 43.99 & \textbf{4.17 $\pm$ 0.32}/ 12.65 & \textbf{22.20 $\pm$ 1.16}/ 14.86\\
			
				EigenAlign & 60.68 $\pm$ 0.29/ 19.37 & 57.44 $\pm$ 0.37/ 19.32 & 66.57 $\pm$ 0.0/ 26.13 & 43.23 $\pm$ 0.0/ 93.60 & 90.51 $\pm$ 0.0/ 2.64 & 98.41 $\pm$ 0.0/ 8.29\\
				FGM & 55.51 $\pm$ 0.0/ 1793.9 & 48.17 $\pm$ 0.0/ 2013.7 & 16.12 $\pm$ 0.0/ 674.94 & \textbf{0.0 $\pm$ 0.0}/ 777.55 & 74.87 $\pm$ 0.05/ 2530.5 & 97.54 $\pm$ 0.01/ 4293.9\\
				LAI-LP & 73.06 $\pm$ 0.23/ 152.47 & 42.00 $\pm$ 0.24/ 154.15 & 49.54 $\pm$ 0.0/ 161.06 & 88.73 $\pm$ 0.1/ 184.153 & 87.98 $\pm$ 0.0/ 33.06 & 96.38 $\pm$ 0.0/ 14.71\\
				PermSync & 10.63 $\pm$ 0.0/ 0.45 & 8.90 $\pm$ 0.0/ 0.46 & 46.93 $\pm$ 0.0/ 0.43 & 79.88 $\pm$ 0.0/ 1.08 & 64.00 $\pm$ 0.0/ 0.22 & 70.00 $\pm$ 0.0/ 0.48\\
			
				RRWM & 60.91 $\pm$ 0.0/ 4.96 & 54.53 $\pm$ 0.0/ 4.83 & 30.99 $\pm$ 0.0/ 8.53 & 1.98 $\pm$ 0.0/ 18.09 & 72.87 $\pm$ 0.01/ 7.98 & 87.04 $\pm$ 0.0/ 21.84\\
				Tensor & 24.37 $\pm$ 0.9/ 93.36 & 15.07 $\pm$ 1.0/ 93.97 & \textbf{1.07 $\pm$ 0.17}/ 107.93 & \textbf{0.0 $\pm$ 0.0}/ 182.07 & 43.24 $\pm$ 2.98/ 40.21 & 32.35 $\pm$ 0.15/ 40.41\\
				IPFP & 65.13 $\pm$ 0.0/ 6.35 & 60.81 $\pm$ 0.0/ 6.28 & 40.90 $\pm$ 0.0/ 8.43 & 3.94 $\pm$ 0.0/ 12.31 & 76.19 $\pm$ 0.0/ 4.65 & 87.74 $\pm$ 0.0/ 8.90\\
			
				PM & 74.63 $\pm$ 0.0/ 7.07 & 71.93 $\pm$ 0.0/ 4.90 & 70.27 $\pm$ 0.0/ 0.94 & 48.82 $\pm$ 0.0/ 1.69 & 83.79 $\pm$ 0.02/ 2.98 & 91.43 $\pm$ 0.24/ 0.44\\
				SMAC & 70.00 $\pm$ 0.0/ 5.75 & 67.36 $\pm$ 0.0/ 5.52 & 50.53 $\pm$ 0.0/ 4.15 & 5.52 $\pm$ 0.0/ 6.90 & 78.56 $\pm$ 0.22/ 1.94 & 87.88 $\pm$ 0.11/ 3.42 \\
				SM & 68.54 $\pm$ 0.0/ 3.34 & 67.18 $\pm$ 0.0/ 3.47 & 65.96 $\pm$ 0.0/ 3.32 & 34.16 $\pm$ 0.0/ 4.93 & 80.27 $\pm$ 0.07/ 1.88 & 88.66 $\pm$ 0.09/ 2.10\\
				GA & 65.06 $\pm$ 0.0/ 4.53 & 64.60 $\pm$ 0.0/ 4.61 & 61.58 $\pm$ 0.0/ 4.08 & 31.62 $\pm$ 0.0/ 5.83 & 77.20 $\pm$ 0.27/ 3.51 & 87.02 $\pm$ 0.16/ 32.28\\
				Munkres & 33.71 $\pm$ 0.0/ 1.52 & 29.99 $\pm$ 0.0/ 1.49 & 51.87 $\pm$ 0.0/ 1.39 & 79.69 $\pm$ 0.0/ 2.45 & 74.00 $\pm$ 0.0/ 0.73 & 92.00 $\pm$ 0.0/ 1.16\\
				
				\hline
		\end{tabular}}
		\hfill{}
		\caption{Car, Motorbike~\cite{cho2013learning}, Butterfly, Magazine~\cite{jiang2011linear}, Building and Books error (\%) for pairwise matchings. Computation time (in seconds) is mentioned after the "/" in the above Table.}
		\label{tb:n2}
	\end{table*}
\end{center}
\vspace{-0.8cm}
\textbf{Rotation:} 
Figure~\ref{fig:first11} shows a $180^{\circ}$ rotated version of the original house frame (Figure~\ref{fig:zero11}). Table \ref{tb:Transformation} shows errors in matching when 
$20$\% of the frames are rotated by both $20^{\circ}$ and $60^{\circ}$, respectively 
and when the same transformations are applied to $40$\% of the frames. As the percentage of transformed frames with greater degree increases, we note a substantial increase in error for other methods in comparison to our method's error increase. 

\textbf{Reflection:} 
The reflected version of a house frame is shown in Figure~\ref{fig:second11}. 
Table \ref{tb:Transformation} shows that affinity-based approaches also performed equally well for reflection of house frame sequences.

\textbf{Scaling:} 
Resizing an image both horizontally and vertically scales the image as is shown in Figure~\ref{fig:third11} . We fixed the scales to $0.5$, $0.75$, $1.25$, and $1.5$ randomly in both the directions in order to transform the images. Our method in Table \ref{tb:Transformation} produces much better matchings than the other methods.

\textbf{Shearing:} 
We randomly apply shearing on house in one of the directions with shear factor $0.5$ (shown in Figure~\ref{fig:fourth11}) and measured the performance shown in Table \ref{tb:Transformation}. In addition to our method, we find that affinity-based algorithms also produce robust matchings.
 
\subsection{Effect of Incomplete and Occluded Landmarks}
To understand the effect of occlusions, we took two real-world datasets, i.e., \emph{Books} and \emph{Building}~\cite{pachauri2013solving} with severe occlusions which are scenes of the same 3D object taken from arbitrary camera angles. These datasets have widely been used in \emph{Structure from Motion (SfM)} problems and are known to be difficult for matching. 
Focusing our attention to the last two columns of Table \ref{tb:n2}, it is evident that our method gives the best results.
\begin{figure}[t!]
	\centering
	\subfigure{%
	
		\includegraphics[trim={0 0 0 0.45cm},clip,width=70mm,height=20mm]{Experiments/Books/BookMatch2_17}}%
	\qquad

	\subfigure{%
	
		\includegraphics[trim={0 0 0 0.45cm},clip,width=70mm,height=20mm]{Experiments/Books/BookMatch5_20}}%
	\caption{Two instances of matchings in Books dataset which is severely occluded. Yellow/green lines show correct/incorrect matches and isolated points show no matches.}
	\label{fig:books}
\end{figure}
Figure~\ref{fig:books} shows the Books dataset where books are placed on a table in various orientations with varying levels of occlusion, along with two sample matchings between different pairs of images. Note that in Figure~\ref{fig:books}, when a corresponding matching clique is not found in the other image, a match isn't \emph{forced} but rather there is no match reported, which doesn't degrade the matching accuracy. Matching as many random cliques, in order of decreasing dimensionality, as possible, manifests itself as an advantage over existing methods, especially when dealing with clutter and/or occlusions.
\begin{figure}[t]
	\centering    
	\includegraphics[width=0.54\linewidth]{Experiments/House_Occlusion_40Plot_Tensor}
	\caption{Error (\%) in matching when varying the number of missing landmarks in	$40\%$ of the images in the frame sequence.}
	\label{fig:occ}
\end{figure}
\textbf{Simulating missing points:}
In order to gain a deeper insight into the behavior of all the matching algorithms, we 
omit $2, 4, 6, 8$, and $10$ ($6.66\%$, $13.33\%$, $20\%$, $26.66\%$, and $33.33\%$) 
points out of total House landmark points (i.e., $30$ points) from $40\%$ (Figure~\ref{fig:occ})
 of frame sequences randomly. In general, all algorithms show an increase in error as more points are removed, but our method has a less gradual increase, while eigenvalue related methods show a rather steep increase in error. Our method is comparable to FGM and RRWM, but the gap in error increases with more missing points. We also observed that FGM incurs the longest runtime for matching in this scenario.

\subsection{Effect of Frame Separation}
Here, we pick two frames from a video for matching and vary the separation in their frame sequence number. The farther apart two frames are the more \emph{pronounced} is the effect we seek between the frame images. For example, as the frame separation increases, \href{http://vasc.ri.cmu.edu/idb/html/motion/hotel/index.html}{\emph{CMU Hotel}} undergoes a more severe 3D rotation, while \emph{Horse-Shear}~\cite{caetano2009learning} undergoes a larger degree of shear. 
\begin{figure}[t]
	\centering
	\parbox{3.6cm}{
		\includegraphics[width=3.6cm]{Experiments/Hotel_FrameSep_Plot_Tensor}}
		\qquad
		\begin{minipage}{3.6cm}
			\includegraphics[width=3.6cm]{Experiments/HorseShear_FrameSep_Plot_Tensor}
		\end{minipage}
		\caption{Error (\%) in matching by various methods with different frame separation level for 
			\emph{CMU Hotel} (on left) and \emph{Horse Shear} (on right).}
			\label{fig:frame-seps}
\end{figure}
We set $p=0.7$ and $k=7$ as nearest neighbors to get the correct matchings.
In Figure~\ref{fig:frame-seps}, in both the left and right plots we notice that most methods show a very sharp rise in error, while our method is quite stable and reports a $0\%$ error. We observe that the naive baseline, LAI-LP also does well and doesn't exhibit steep changes in error with larger frame separation.

\textbf{Experimental Summary:}
In general, we find that the \emph{affinity matrix} based methods like FGM and RRWM are more robust to affine transformations than other competing algorithms. Our method performs the best as the weight vectors in our algorithm effectively capture even the \emph{higher-order geometric properties} of the neighborhood and nearly preserves them under affine transformations. The naive baseline, i.e., LAI-LP, does not perform as well because it also has a feature-based component like SIFT which is known to fail on some affine transformations.

\section{Conclusion}
To the best of our knowledge, we have presented the first approach towards partial higher-order matching by initially capturing higher-order structure as \emph{random clique complexes} and then proposing a corresponding matching algorithm. 
From a theoretical point of view, we studied matching as a QAP on random clique adjacency matrices that represented the $k$-skeleton of our random clique complexes and gave bounds on the concentration inequality of the spread of its eigenvalues. We also improved bounds on the largest eigenvalue of the Lawler QAP formulation, used by affinity-matrix based approaches. We discussed the robustness of such approaches to missing points and also showed the sharpness of our result. Furthermore, inspired by Finke et. al.~\cite{martello1987} we studied the asymptotic behavior of our higher dimensional clique adjacency matrices. A more detailed investigation of the distribution of eigenvalue gaps for such random matrices with Poisson distributed entries is left for future work. 

From an empirical perspective, we compared the matching accuracies of diverse algorithms on both synthetic and real-world datasets that were known to have severe occlusions and distortions, thus posing a daunting challenge to matching algorithms. We argue that our experiments show strong evidence that our approach outperforms all the state-of-the-art matching methods on a diverse range of datasets.  

\clearpage
\section*{Acknowledgements}
We thank our colleagues from the Mathematics Dept. at IIT-H (Sukumar Daniel, Narasimha Kumar, and Bhakti B. Manna) for their insight and expertise.
We would also like to thank all the reviewers for their feedback and suggestions. We are grateful to the authors of ~\cite{ZhouD16,zhou2013deformable,cho2010reweighted,feizi2016spectral,pachauri2013solving,li2013object,duchenne2011tensor} for providing their source codes and datasets.

