Solving Partial Assignment Problems using Random Clique Complexes

TL;DR

This paper introduces a novel approach to partial assignment problems by leveraging random clique complexes, which capture higher-order structures, and demonstrates superior accuracy and robustness over existing methods through theoretical analysis and experiments.

Contribution

It proposes a new formulation of partial assignment problems using random clique complexes and provides theoretical and empirical validation of its effectiveness.

Findings

Outperforms existing matching algorithms significantly.

Effective on datasets with severe occlusions and distortions.

Provides theoretical analysis of runtime and asymptotic behavior.

Abstract

We present an alternate formulation of the partial assignment problem as matching random clique complexes, that are higher-order analogues of random graphs, designed to provide a set of invariants that better detect higher-order structure. The proposed method creates random clique adjacency matrices for each k-skeleton of the random clique complexes and matches them, taking into account each point as the affine combination of its geometric neighbourhood. We justify our solution theoretically, by analyzing the runtime and storage complexity of our algorithm along with the asymptotic behaviour of the quadratic assignment problem (QAP) that is associated with the underlying random clique adjacency matrices. Experiments on both synthetic and real-world datasets, containing severe occlusions and distortions, provide insight into the accuracy, efficiency, and robustness of our approach. We outperform diverse matching algorithms by a significant margin.

Tables

Table 1. Table 1: Error (%) of transformation on CMU House: inserted 20 % percent 20 20\% and 40 % percent 40 40\% impurity in CMU House frame sequence randomly for rotation ( 20 ∘ superscript 20 20^{\circ} , 60 ∘ superscript 60 60^{\circ} ), reflection, scaling and shear. Minimum error (%) is shown in bold. Matching is computed for 111 111 111 frames from the 1 s ​ t superscript 1 𝑠 𝑡 1^{st} frame to the other 110 110 110 frames. Our method shows best performance among all the methods.

Algorithms	${20}^{\circ }$ Rotation		${60}^{\circ }$ Rotation		Reflection		Scaling		Shear
	20%	40%	20%	40%	20%	40%	20%	40%	20%	40%
OurMethod	0.01 $\pm$ 0.0	0.01 $\pm$ 0.0	0.05 $\pm$ 0.0	0.03 $\pm$ 0.0	0.0 $\pm$ 0.0	0.0 $\pm$ 0.0	0.1 $\pm$ 0.0	0.2 $\pm$ 0.1	0.0 $\pm$ 0.0	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	0.0 $\pm$ 0.0	0.0 $\pm$ 0.0	2.20 $\pm$ 1.0	3.4 $\pm$ 0.5	0.0 $\pm$ 0.0	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	0.0 $\pm$ 0.0	0.0 $\pm$ 0.0	5.0 $\pm$ 2.0	10.0 $\pm$ 1.5	0.0 $\pm$ 0.0	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	0.0 $\pm$ 0.0	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

Table 2. Table 2: Car, Motorbike  (Cho et al., 2013 ) , Butterfly, Magazine  (Jiang et al., 2011 ) , Building and Books error (%) for pairwise matchings. Computation time (in seconds) is mentioned after the "/" in the above Table.

Methods	Car	Bike	Butterfly	Magazine	Building	Book
OurMethod	4.14 $\pm$ 2.45/ 7.13	3.15 $\pm$ 0.32/ 6.96	3.89 $\pm$ 0.23/ 14.76	0.48 $\pm$ 0.02/ 43.99	4.17 $\pm$ 0.32/ 12.65	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	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	1.07 $\pm$ 0.17/ 107.93	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

Table 3. Table 3: Neighbourhood of 3 , 2 , 1 3 2 1 3,2,1 -cliques of graphs G 1 subscript 𝐺 1 G_{1} and G 2 subscript 𝐺 2 G_{2} shown in Figure  6 .

Cliques	Graph 1 ($G_1$)		Graph 2 ($G_2$)
	Clique	Neighbours	Clique	Neighbours
3-Cliques	{A,B,D}	{A}, {B}, {D}, {A,B}, {A,D}, {B,D}, {B,D,E}	{A,C,D}	{A}, {C}, {D}, {A,C}, {A,D}, {C,D}, {B,D,E}
	{B,D,E}	{B}, {D}, {E}, {B,D}, {B,E}, {D,E}, {A,B,D}	{B,D,E}	{B}, {D}, {E}, {B,D}, {B,E}, {D,E}, {A,C,D}
2-Cliques	{A,B}	{A}, {B}, {A,B,D}	{A,C}	{A}, {C}, {A,C,D}
	{A,C}	{A}, {C}	{A,D}	{A}, {D}, {A,C,D}
	{A,D}	{A}, {D}, {A,B,D}	{B,D}	{B}, {D}, {B,D,E}
	{B,D}	{B}, {D}, {A,B,D}, {B,D,E}	{B,E}	{B}, {E}, {B,D,E}
	{B,E}	{B}, {E}, {B,D,E}	{C,D}	{C}, {D}, {A,C,D}
	{D,E}	{D}, {E}, {B,D,E}	{D,E}	{D}, {E}, {B,D,E}
	{D,F}	{D}, {F}	{D,F}	{D}, {F}
1-Cliques	{A}	{A,B}, {A,C}, {A,D}, {A,B,D}	{A}	{A,C}, {A,D}, {A,C,D}
	{B}	{A,B}, {B,D}, {B,E}, {A,B,D}, {B,D,E}	{B}	{B,D}, {B,E}, {B,D,E}
	{C}	{A,C}	{C}	{A,C}, {C,D}, {A,C,D}
	{D}	{A,D}, {B,D}, {D,E}, {D,F}, {A,B,D},	{D}	{A,D}, {B,D}, {C,D}, {D,E}, {D,F}, {A,C,D},
		{B,D,E}		{B,D,E}
	{E}	{B,E}, {D,E}, {B,D,E}	{E}	{B,E}, {D,E}, {B,D,E}
	{F}	{D,F}	{F}	{D,F}

Table 4. Table 4: Matchings of 3 , 2 , 1 3 2 1 3,2,1 -cliques of graphs G 1 subscript 𝐺 1 G_{1} and G 2 subscript 𝐺 2 G_{2} shown in Figure  6 .

Cliques	Graph 1 ($G_1$)	Graph 2 ($G_2$)
3-Cliques	{B,D,E} $\to$ 1	{B,D,E} $\to$ 1
2-Cliques	{A,C} $\to$ 1	{A,C} $\to$ 1
	{A,D} $\to$ 2	{A,D} $\to$ 2
	{B,D} $\to$ 5	{B,D} $\to$ 5
	{B,E} $\to$ 6	{B,E} $\to$ 6
	{D,E} $\to$ 7	{D,E} $\to$ 7
	{D,F} $\to$ 8	{D,F} $\to$ 8
1-Cliques	{A}, {B}, {C}	{A}, {B}, {C}
	{D}, {E}, {F}	{D}, {E}, {F}

Table 5. Table 5: Datasets used, where N 𝑁 N is the number of samples and n 𝑛 n is the dimensionality of each sample.

Groups	Dataset	$N\times n$
Video Frames	CMU House	$111\times 30$
	CMU Hotel	$101\times 30$
Affine	Horse-Rot (Caetano et al., 2009)	$200\times 35$
	Horse-Shear (Caetano et al., 2009)	$200\times 35$
Occluded	Books (Pachauri et al., 2013)	$20\times 34$
	Building (Pachauri et al., 2013)	$16\times 28$
Non-Affine	Magazine (Jiang et al., 2011)	$30\times 30$
	Butterfly (Jiang et al., 2011)	$30\times 19$
Object Matching	Car (Cho et al., 2013)	$40\times 10$
	Bike (Cho et al., 2013)	$40\times 10$

Equations

$$\displaystyle\mathcal{X}^{(k)}(G):=\left(\mathcal{X}^{(k-1)}(G)\cup\coprod_{\sigma:\textit{dim }\sigma=k}\sigma^{(k)}\right)\Bigg{/}\sim$$

$$\underset{X_0,\dots ,X_h}{\mathrm{argmin}}$$

$$\forall k\le h,\text{mathds}{1}^T{X}_k=\text{mathds}1,X_k^T\text{mathds}1=\text{mathds}1$$

$$\text{minimize }\mathrm{Tr}(AXB{X}^T)\text{ s.t. }X\in {\Pi }_X$$

$${\delta }_i(A)\le 2\left|\left|\left|A\right|\right|\right|$$

$$\mathbb{P}(|{\mathcal{S}}_n(A)-\mu |\ge ϵ\sigma )\le K\max \left(e^{(-p{ϵ}^2)},e^{(-pϵ\sigma /2|||A|||)}\right)$$

$$\mathbb{P}\left\{\frac{\underset{\pi \in \Pi }{\max }\sum _{v{v}^{\prime }}{\mathcal{C}}_{v{v}^{\prime }}}{\underset{\pi \in \Pi }{\min }\sum _{v{v}^{\prime }}{\mathcal{C}}_{v{v}^{\prime }}}\le 1+ϵ\right\}$$

$$\ge 1-2|\Pi |\exp \left(-2|S_{\pi }|{\left(\frac{{ϵ}^{\prime }\sqrt{{\lambda }_v}}{{ϵ}^{\prime }+2{\lambda }_v}\right)}^2\right)=:\psi (n,ϵ)$$

$$\text{mathds}{1}_A=\{\begin{array}{ll}1& \text{if A is a k-clique in G(n,p)}\\ 0& \text{otherwise}\end{array}$$

$$\mathbb{E}(X_n(k))=\mathbb{E}(\sum _{|A|=k}\text{mathds}{1}_A)=\sum _{|A|=k}\mathbb{E}(\text{mathds}{1}_A)=\left(\genfrac{}{}{0ex}{}{n}{k}\right)p^{\left(\genfrac{}{}{0ex}{}{k}{2}\right)}$$

$$\sum _{v,v^{\prime }}{b}_{v{v}^{\prime }}{a}_{\pi (v)\pi (v^{\prime })}$$

$$\pi \in \Pi$$

$$S_{\pi }=\{(\pi (v),\pi (v^{\prime })\mid v<v^{\prime },u,v=1,\dots ,n\}$$

$$\mathbb{P}\left\{\frac{1}{|S_{\pi }|}\left|\sum _{v,v^{\prime }\in \pi }({\mathcal{C}}_{v{v}^{\prime }}-{\lambda }_e)\right|\ge {ϵ}^{\prime }\right\}$$

$$D:=\sqrt{\sum _{k=1}^n\mathbb{V}(X_k)}$$

$$\mathbb{P}\left\{\left|\sum _{k=1}^n(X_k-\mathbb{E}(X_k))\right|\ge \mu D\right\}\le 2\exp \left(-\frac{{\mu }^2}{2(1+\mu /2D{)}^2}\right)$$

$$\mathbb{P}\left\{\frac{1}{|S_{\pi }|}\left|\sum _{v,v^{\prime }\in \pi }({\mathcal{C}}_{v{v}^{\prime }}-{\lambda }_e)\right|\ge {ϵ}^{\prime }\right\}$$

$$\le$$

$$D=\sqrt{\sum _{v,v^{\prime }\in \pi }{\lambda }_v}=\sqrt{|S_{\pi }|{\lambda }_v}$$

$$|\Pi |\mathbb{P}\left\{\left(\frac{\sqrt{{\lambda }_v}}{\sqrt{|S_{\pi }|}}\right)\left(\frac{1}{\sqrt{|S_{\pi }|{\lambda }_v}}\right)\left|\sum _{v,v^{\prime }\in \pi }({\mathcal{C}}_{v{v}^{\prime }}-{\lambda }_e)\right|\ge {ϵ}^{\prime }\right\}$$

$$=$$

$$|\Pi |\mathbb{P}\left\{\left|\sum _{v,v^{\prime }\in \pi }({\mathcal{C}}_{v{v}^{\prime }}-{\lambda }_e)\right|\ge \left(\frac{{ϵ}^{\prime }\sqrt{|S_{\pi }|}}{\sqrt{{\lambda }_v}}\right)\left(\sqrt{|S_{\pi }|{\lambda }_v}\right)\right\}$$

$$\le$$

$$=$$

$$\mathbb{P}\left\{\frac{1}{|S_{\pi }|}\left|\sum _{v,v^{\prime }\in \pi }({\mathcal{C}}_{v{v}^{\prime }}-{\lambda }_e)\right|\le {ϵ}^{\prime }\right\}\ge 1-$$

$$2|\Pi |\exp \left(-2|S_{\pi }|{\left(\frac{{ϵ}^{\prime }\sqrt{{\lambda }_v}}{{ϵ}^{\prime }+2{\lambda }_v}\right)}^2\right)$$

$$|S_{\pi }|({\lambda }_e-{ϵ}^{\prime })\le \sum _{v,v^{\prime }\in \pi }{\mathcal{C}}_{v{v}^{\prime }}\le |S_{\pi }|({\lambda }_e+{ϵ}^{\prime })$$

$$\frac{\underset{\pi \in \Pi }{\max }\sum _{v{v}^{\prime }}{\mathcal{C}}_{v{v}^{\prime }}}{\underset{\pi \in \Pi }{\min }\sum _{v{v}^{\prime }}{\mathcal{C}}_{v{v}^{\prime }}}\le \frac{|S_{\pi }|({\lambda }_e+{ϵ}^{\prime })}{|S_{\pi }|({\lambda }_e-{ϵ}^{\prime })}\le 1+ϵ$$

$$\sum _{i:x_i\ne y_i}|{\alpha }_i|\le t{\left(\sum _{i=1}^m{\alpha }_i^2\right)}^{1/2}$$

$$\mathbb{P}[\mathcal{A}]\mathbb{P}[{\bar{\mathcal{A}}}_t]\le e^{-t^2/4}.$$

$$\mathbb{P}[|{\lambda }_1(A)-\mathcal{M}|\ge t]\le 4e^{-t^2/8},$$

$$\mathcal{R}(A,x)=\frac{x^{T}Ax}{x^{T}x}$$

$${\lambda }_1(A)=\mathcal{R}(A,v)=\frac{v^{T}Av}{v^{T}v}$$

$${\lambda }_1(A)=\mathcal{R}(A,v)=v^{T}Av=\underset{\text{off-diagonal}}{\underbrace{\sum _{1\le i<j\le m}(v_i^T{v}_j+v_j^T{v}_i)a_{ij}}}+$$

Full text

Machine Learning, ICML

1 Introduction

The 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 (Singer & Shkolnisky, 2011), protein structure comparisons in bioinformatics (Zaslavskiy et al., 2009), and computer vision (Conte et al., 2004). Computer vision especially boasts a broad range of applications that include object matching, image registration (Shen & Davatzikos, 2002), stereo matching (Goesele et al., 2007), shape matching (Petterson et al., 2009; Berg et al., 2005), structure from motion (SfM) (Szeliski, 2010), and object detection (Jiang et al., 2011), 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 linear assignment problem (LP) versus ones that solve a quadratic assignment problem (QAP) via 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 relational structure is best preserved.

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.

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 (Chertok & Keller, 2010).

Motivated by the aforementioned observations and inspired by Kahle (Kahle, 2006)’s work on combinatorial topological models like the random clique complex, we focus our attention to matching higher-order components between two sets of points in the setting of some points 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 1 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 clique complex is the simplicial complex with all complete subgraphs (i.e., cliques) as its 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 (Giusti et al., 2015). 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.

Contributions: We present the following contributions.

To the best of our knowledge, our proposed approach is a first attempt to formulate higher-order matching between two sets of points, given partial or incomplete information, as a matching between two random clique complexes. We also propose an efficient matching algorithm and study both its time and storage complexity. 2. 2.

(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 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 clique percolation regime (Bollobás & Riordan, 2009), where the entries follow a Poisson distribution. 3. 3.

Finally, we present a comprehensive empirical study that compares our method’s matching accuracy to that of a diverse set of matching approaches (Zhou & De la Torre, 2016; Zhou & De la Torre, 2013; Cho et al., 2010; Feizi et al., 2016; Leordeanu & Hebert, 2005; Cour et al., 2007; Pachauri et al., 2013; Gold & Rangarajan, 1996; Kuhn, 1955; Leordeanu et al., 2009; Zass & Shashua, 2008; Li et al., 2013; Duchenne et al., 2011). 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.

2 Matching Random Clique Complexes

We consider the problem of capturing higher-order feature groups among landmark points in an image by representing them as a 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.

2.1 Structure of a Random Clique Complex

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^{\prime}\}\in{V\choose 2}$, are i.i.d Bernoulli($p$) distributed. Recall, that a $k$-clique in $G(n,p)$ is a 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 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 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

[TABLE]

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)$.

$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ás et. al. (Bollobás & Riordan, 2009), 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 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.

2.2 Problem 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^{\prime})$, for all $i\leq h$, subject to assignment constraints. This can be formulated as a 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^{\prime})=\{G^{\prime(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^{\prime(k,l)}$. The combinatorial matching requires the optimal set of permutation matrices that best align $\mathcal{X}(G)$ and $\mathcal{X}(G^{\prime})$. More formally, this can be expressed as the following constrained optimization problem

[TABLE]

2.3 Our Algorithm

At a high level, our goal is to minimize $\lVert\mathcal{X}(G)-\mathcal{X}(G^{\prime})\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 (1) aims to find a set of permutation matrices $\{X_{1},\dots,X_{h}\}$ that minimizes the overall number of misalignments between equi-dimensional faces of $\mathcal{X}(G)$ and $\mathcal{X}(G^{\prime})$, 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^{\prime})$.

Algorithm 1 presents our method to solve the combinatorial optimization problem (Equation (1)). In decreasing order of clique dimensionality, for a fixed dimension $k$ and given the adjacency matrices $G^{(k,l)}$ and $G^{\prime(k,l)}$ for $k$-skeletons $\mathcal{X}^{(k)}(G)$ and $\mathcal{X}^{(k)}(G^{\prime})$, respectively. In every iteration, our objective is to solve $\operatornamewithlimits{argmin}_{X_{k}}\lVert G^{(k,l)}X_{k}-X_{k}G^{\prime(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 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 geometric properties of the neighborhood of every clique. We achieve this by characterizing the $i$-th clique’s barycenter $c_{i}^{(k)}$ as an 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 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^{\prime}$ (Step $13$). Finally, the Kuhn-Munkres (Kuhn, 1955) 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.

2.4 Complexity Analysis

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 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 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ás (Bollobás & Riordan, 2009), where they identify a threshold probability for 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 Poisson distribution. Exceeding this threshold results in formation of 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 $knp^{k}$ expected relocations, which in turn estimates $nnz(G^{(k,l)})$.

Note that for every iteration in Algorithm 1, the dominating cost is that of running the Kuhn-Munkres matching algorithm in Step $14$, which has a 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})$.

3 Theoretical Analysis of QAP

In this section, we present three related results in the context of matching random matrices, namely: (i) concentration inequality of eigenvalue bounds on the QAP trace formulation for random symmetric matrices, (ii) tighter concentration inequality of eigenvalue bounds on the Lawler QAP formulation on random symmetric matrices in the context of works that use 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$.

3.1 Eigenvalue Bounds of Trace QAP Formulation on Random Matrices

Let $A=(a_{vv^{\prime}})$, $B=(b_{vv^{\prime}})\in\mathds{R}^{n\times n}$ be random real-symmetric matrices. Let $X=(x_{ij})\in\mathds{R}^{n\times n}$ be a permutation matrix. Then, the trace formulation of a QAP is given by

[TABLE]

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)$111The 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}=\operatorname{diag}(\lambda_{1}(A),\dots,\lambda_{n}(A))$ and $\Lambda_{B}=\operatorname{diag}(\lambda_{1}(B),\dots,\lambda_{n}(B))$ with their corresponding orthogonal eigenvector matrices $Q_{A}$ and $Q_{B}$. Finke et. al. (Martello et al., 1987) gave the following eigenvalue bounds.

Theorem 1**.**

Let $A$ and $B$ be symmetric matrices. Then for all $X\in X_{\Pi}$,

$\operatorname{Tr}(AXBX^{T})=\lambda(A)^{T}Q^{(X)}\lambda(B)$*, *

where $Q^{(X)}=\langle Q_{A}^{(i)},XQ_{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; 2. 2.

$\mathcal{L}\leq\operatorname{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)$* *

It was also noticed by Finke (Martello et al., 1987) that these bounds can further be tightened by reducing the 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. (Martello et al., 1987) suggested a reduction method to further sharpen the bound by replacing matrices $A$ and $B$ by smaller spread symmetric matrices $\widetilde{A}$ and $\widetilde{B}$. The reductions are achieved as $\widetilde{A}=A-M_{A}-M_{A}^{T}-\mathfrak{D}_{A}$ and $\widetilde{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}(\widetilde{A})$ and $\mathfrak{S}(\widetilde{B})$.

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 $\widetilde{A}\in\mathds{R}^{n\times n}$, with eigenvalues $\lambda_{1}(\widetilde{A})\leq\dots\leq\lambda_{n}(\widetilde{A})$, we define the gap (spacing) between its consecutive eigenvalues as $\delta_{i}(\widetilde{A}):=|\lambda_{i+1}(\widetilde{A})-\lambda_{i}(\widetilde{A})|$ for $1\leq i\leq n-1$. Then, the spread $\mathfrak{S}(\widetilde{A})$ for the reduced matrix $\widetilde{A}$ can be redefined as: $\mathfrak{S}(\widetilde{A})=\sum_{i=1}^{n}\delta_{i}(\widetilde{A})$

We begin by upper bounding $\delta_{i}(\widetilde{A})$ using the following lemma 1 (proof in supplementary notes).

Lemma 1**.**

Let ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|.\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ 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}$,

[TABLE]

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{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|A\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$, as shown in Lemma 1. 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{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|A\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ a.s, with mean $\mu_{i}(A)$ and variance $\sigma_{i}^{2}(A)$. Then, using Chernoff bounds, for any $\epsilon>0$, we have

[TABLE]

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.

3.2 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 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 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. (Leordeanu & Hebert, 2005) show via a spectral relaxation that Lawler’s WAP reduces to solving $w^{*}=\mathop{\mathrm{argmax}}\nolimits_{w}\frac{w^{T}Aw}{w^{T}w}$, $w\in\mathds{R}^{n^{2}}$. This is solved by finding the leading eigenvalue $\lambda_{1}(\mathscr{A})$.

As illustrated in (Alon et al., 2002), we also use Talagrand’s concentration inequality (Talagrand, 1995). We provide a tighter bound in the case of our affinity matrix using Rayleigh’s quotient.

Theorem 2**.**

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})$.∎

Discussion: We further investigate the robustness of affinity-matrix based graph matching solutions when dealing with missing or incomplete data. We show the sharpness of our result in Theorem 5 on the affinity matrix, similar to (Alon et al., 2002), who analyze their results using 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 (Alon et al., 2002), 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 average degree of $G$ is given by $\bar{d}=\sum_{v\in V}d_{G}(v)/|V|$ and its 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 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 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 tight up to a constant factor in the exponent. Our experimental results on Factorized Graph Matching (FGM) by Zhou et. al. (Zhou & De la Torre, 2016) and Re-weighted Random Walk Matching (RRWM) (Cho et al., 2010) also support the finding that affinity matrix based matching solutions are more robust to missing edges due to occlusions in data.

3.3 Asymptotic Analysis of Higher-Order Clique Assignment

Following along the same lines as Finke et. al. (Martello et al., 1987), we study the asymptotic behavior of the worst to most optimal ratio and present it as the following theorem.

Theorem 3**.**

Given random clique adjacency matrices $A^{k,l}_{n},B^{k,l}_{n}\sim Pois(\lambda)$ and their associated cost matrix $\mathcal{C}_{vv^{\prime}}\sim Pois(\lambda)$. We denote by $\lambda_{e}:=\mathbb{E}(\mathcal{C}_{vv^{\prime}})$ and $\lambda_{v}:=\mathbb{V}(\mathcal{C}_{vv^{\prime}})$ 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

[TABLE]

where, $|\Pi|=n!$, $|S_{\pi}|={{n+1}\choose 2}$, $\lim_{n\to\infty}\psi(n,\epsilon)=1$∎

4 Experiments

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) affinity-matrix: FGM (Zhou & De la Torre, 2016; Zhou & De la Torre, 2013)222FGM, RRWM (Cho et al., 2010), (ii) Eigenvalues: EigenAlign (Feizi et al., 2016)333EigenAlign, SM (Leordeanu & Hebert, 2005), SMAC (Cour et al., 2007), PermSync (Pachauri et al., 2013)444PermSync, (iii) LP relaxation: GA (Gold & Rangarajan, 1996), Kuhn-Munkres (Kuhn, 1955), (iv) Integer QAP: IPFP (Leordeanu et al., 2009), (v) Probabilistic matching: PM (Zass & Shashua, 2008), (vi) Higher-order matching given complete data: Tensor (Duchenne et al., 2011)555Tensor, and (vii) Geometric and Feature matching: LAI-LP (Li et al., 2013)666This algorithm serves as our naive baseline as it directly uses neighborhood properties of the underlying graph (LAI-LP).. Our code777Our Method is publicly available.

4.1 Effect of Affine Transformations

Simulated Dataset: We perform affine transformations on 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 2 shows examples of affine transformations on house frame sequences. Table 1 shows the comparative error in matching for all the algorithms. We now describe each affine transformation as performance metrics in our experiments.

Rotation: Figure 8 shows a $180^{\circ}$ rotated version of the original house frame (Figure 8). Table 1 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.

Reflection: The reflected version of a house frame is shown in Figure 8. Table 1 shows that affinity-based approaches also performed equally well for reflection of house frame sequences.

Scaling: Resizing an image both horizontally and vertically scales the image as is shown in Figure 8 . 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 1 produces much better matchings than the other methods.

Shearing: We randomly apply shearing on house in one of the directions with shear factor $0.5$ (shown in Figure 8) and measured the performance shown in Table 1. In addition to our method, we find that affinity-based algorithms also produce robust matchings.

4.2 Effect of Incomplete and Occluded Landmarks

To understand the effect of occlusions, we took two real-world datasets, i.e., Books and Building (Pachauri et al., 2013) with severe occlusions which are scenes of the same 3D object taken from arbitrary camera angles. These datasets have widely been used in Structure from Motion (SfM) problems and are known to be difficult for matching. Focusing our attention to the last two columns of Table 2, it is evident that our method gives the best results.

Figure 3 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 3, when a corresponding matching clique is not found in the other image, a match isn’t 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.

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 4) 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.

4.3 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 pronounced is the effect we seek between the frame images. For example, as the frame separation increases, CMU Hotel undergoes a more severe 3D rotation, while Horse-Shear (Caetano et al., 2009) undergoes a larger degree of shear.

We set $p=0.7$ and $k=7$ as nearest neighbors to get the correct matchings. In Figure 5, 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.

Experimental Summary: In general, we find that the 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 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.

5 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 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. (Martello et al., 1987) 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.

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  (Zhou & De la Torre, 2016; Zhou & De la Torre, 2013; Cho et al., 2010; Feizi et al., 2016; Pachauri et al., 2013; Li et al., 2013; Duchenne et al., 2011) for providing their source codes and datasets.

Appendix A Proofs

A.1 Upper Bound to Clique Size in a Random Graph

Let $G(n,p)$ denote the Erdős-Rényi random graph on $n$ vertices, i.e., $G(n,p)=\{G_{ij}|1\leq i<j\leq n\}$, where $G_{ij}\sim Ber(p)$ are i.i.d Bernoulli random variables. We denote the number of $k$-cliques in the realization of $G(n,p)$ as $X_{n}(k)$. By definition, a $k$-clique in a graph $G$ is a subset $A$ of $k$ vertices, which induce a complete subgraph of $G$. Additionally, no other vertex in $G$ can be joined by edges to all vertices of $A$. Therefore, we can represent $X_{n}(k)$ as a sum of indicator random variables $\mathds{1}_{A}$, where

[TABLE]

It is clear that $X_{n}(k)=\sum_{|A|=k}\mathds{1}_{A}$. Hence, we get

[TABLE]

Using Stirling’s formula, we upper bound $X_{n}(k)$ as $\left(\frac{en}{k}\right)^{k}$, where $e$ is the Euler’s number.

A.2 Quadratic Assignment Problem

We begin by defining the general quadratic assignment problem (QAP) using the Koopman-Beckmann version. Let $A=(a_{vv^{\prime}})$, $B=(b_{vv^{\prime}})\in\mathds{R}^{n\times n}$. Let $\Pi$ denote the set of all possible bijections (permutations) $\pi:N\rightarrow N$, where $N=\{1,2,\dots n\}$. We define the QAP as:

[TABLE]

For now on, for ease of notation, we denote the cost function $b_{vv^{\prime}}a_{\pi(v)\pi(v^{\prime})}$ as $\mathcal{C}_{vv^{\prime}}$.

A.3 Asymptotic Analysis of Higher-order Clique Assignment (Proof of Theorem $3$)

Given that the QAP is a combinatorial optimization problem, in the case of random symmetric matrices, the subset of feasible solutions $S_{\pi}$ is of the form:

[TABLE]

where, $|S_{\pi}|={n+1\choose 2}$ and $|\Pi|=n!$.

Recall that our cost function $\mathcal{C}_{vv^{\prime}}$ has expectation $\lambda_{e}$ and variance $\lambda_{v}$. For notational convenience, we set $\epsilon^{\prime}=\lambda_{v}-\epsilon$ . Then, there exists a bijection $\pi\in\Pi$, for which the following holds by the definition of variance

[TABLE]

To proceed further with our proof, we make use of the following lemma by Renyi et. al. (Rényi, 1970).

Lemma 2**.**

Let $X_{1},\dots,X_{n}$ be independent random variables with $|X_{k}-\mathbb{E}(X_{k})|\leq 1$, $k=1,\dots,n$. Denote

[TABLE]

and let $\mu$ be a positive real number with $\mu\leq D$. Then

[TABLE]

∎

In order to apply Lemma 2, we change the form of the inequality as follows:

[TABLE]

Before applying Lemma 2, we compute $D$ as,

[TABLE]

We can rewrite (5) as

[TABLE]

Now, we make use of Lemma 2 and get

[TABLE]

Equation 3 can now be written as

[TABLE]

It can easily be verified that the expression in the R.H.S. of the above inequality tends to $1$ as $n\to\infty$.

We know that for the expression $\left|\sum_{v,v^{\prime}\in\pi}(\mathcal{C}_{vv^{\prime}}-\lambda_{e})\right|\leq\epsilon^{\prime}|S_{\pi}|$, the following bounds hold.

[TABLE]

It follows that

[TABLE]

This completes the proof.∎

A.4 Eigenvalue Bounds on Lawler’s QAP Formulation on Random Matrices (Proof of Theorem 2)

As illustrated in (Alon et al., 2002), we will make use of Talagrand’s concentration inequality. We provide a tighter bound in the case of our affinity matrix using the Rayleigh’s quotient.

Theorem 4**.**

(Talagrand, 1995)*

Let $\Omega=\prod_{i=1}^{m}\Omega_{i}$ be a product space of probability spaces. Let $\mathcal{A}$ and $\mathcal{A}_{t}$ be subsets of $\Omega$ and if for each $y=(y_{1},\dots,y_{m})\in\mathcal{A}_{t}$, there exists a real vector $\alpha=(\alpha_{1},\dots,\alpha_{m})$, such that for every $x=(x_{1},\dots,x_{k})\in\mathcal{A}$, the following inequality holds*

[TABLE]

Then,

[TABLE]

Here, $\mathcal{A}_{t}$ denotes the set with Talagrand distance at most $t$ from $\mathcal{A}$ and $\bar{\mathcal{A}}_{t}$ denotes the complement of set $\mathcal{A}_{t}$.∎

Theorem 5**.**

For a real symmetric matrix $A=(a_{ij})\in R^{m\times m}$ and for positive constant $t$,

[TABLE]

where $\mathcal{M}$ is the median of $\lambda_{1}(A)$.

Proof.

888Our proof technique follows the technique outlined in (Alon et al., 2002)

Given a real symmetric matrix $A=(a_{ij})\in R^{m\times m}$ and a non-zero vetor $x$, the Rayleigh Quotient $\mathcal{R}(A,x)$ is defined as

[TABLE]

Given the eigenvalues of $A$ in decreasing order as $\lambda_{1}(A)\geq\dots\geq\lambda_{m}(A)$, we know that $\mathcal{R}(A,x)\in[\lambda_{m}(A),\lambda_{1}(a)]$. It is well known that $\mathcal{R}(A,x)$ attains its maximum value at $\lambda_{1}(A)$ when $x=v$, where $v$ is the eigenvector corresponding to $\lambda_{1}(A)$. Therefore, we have

[TABLE]

In our proof, we omit the constant factor $v^{T}v$ and normalize the eigenvector $v$, hence $\left\lVert v\right\rVert=1$.

Consider the product space $\Omega$ of entries $a_{ij}$, $1\leq i\leq j\leq m$. Let $t,\mathcal{M}$ be real numbers, where $t>0$ and $\mathcal{M}$ is the median of $\lambda_{1}(A)$. Let $\mathcal{A}$ be the set of matrices $A=(a_{ij})\in\Omega$, for which $\lambda_{1}(A)\leq\mathcal{M}$. By definition, $\mathbb{P}[\mathcal{A}]\geq 1/2$. Additionally, let $\mathcal{B}$ be the set of matrices $B=(b_{ij})\in\Omega$, for which $\lambda_{1}(B)\geq\mathcal{M}+t$. Using Rayleigh’s equation (6) for $\lambda_{1}(A)$, we rewrite it as a summation of diagonal and off-diagonal terms

[TABLE]

and

[TABLE]

In order to apply Talagrand’s inequality (Theorem 4), we set a real vector $\alpha=(\alpha_{ij})_{1\leq i\leq j\leq m}$ as follows: For off-diagonal $(1\leq i<j\leq m)$ terms, we set

[TABLE]

For diagonal $(1\leq i\leq m)$ terms, we set

[TABLE]

We proceed by first proving two claims that will be used in this proof.

Claim 1**.**

[TABLE]

Proof.

By definition,

[TABLE]

This completes the proof. ∎

Claim 2**.**

For every $A\in\mathcal{A}$,

[TABLE]

Proof.

Recall that for matrix $A\in\mathcal{A}$, $v$ is the eigenvector with unit-norm corresponding to $\lambda_{1}(A)$. We know that,

[TABLE]

while,

[TABLE]

We observe that the entries in affinity matrices $A$ and $B$, are affinity scores in interval $[0,1]$. Therefore, we have $|b_{ij}-a_{ij}|\leq 1$, for all $1\leq i,j\leq m$. For ease of notation, let us denote by $P$, the set of ordered pairs $ij$ with $1\i,j\leq m$ where $a_{ij}\neq b_{ij}$. Then,

[TABLE]

This completes the proof. ∎

By the above two claims, we get the following form:

[TABLE]

Applying Talagrand’s inequality, we get

[TABLE]

Since $\mathcal{M}$ is the median of $\lambda_{1}(A)$, by definition $\mathbb{P}[\lambda_{1}(A)\leq\mathcal{M}]\geq 1/2$, then

[TABLE]

Accordingly, we also have that,

[TABLE]

Combining results (7) and (8), we have

[TABLE]

This completes the proof. ∎

A.5 Proof of Lemma 1

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 $Ax_{i}=\lambda_{i}x_{i}$, we have

[TABLE]

It follows,

[TABLE]

As ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|X\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ is non-negative, we get $|\lambda_{i}|\leq{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|A\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$. Thus, every eigenvalue of $A$ is upper bounded by the matrix norm ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|A\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$. Applying the triangle inequality, we get that $\delta_{i}(A)=|\lambda_{i+1}(A)-\lambda_{i}(A)|\leq 2{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|A\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$, which completes the proof.∎

Appendix B Example

We explained our method with the help of an example shown in Figure 6, Table 3 and Table 4 for a better understanding. We consider two random graphs $G_{1}$ and $G_{2}$ with $6$ vertices each in Figure 6 for which we perform higher-order matching from $3$-cliques 6 to $1$-cliques 6. For a higher-order matching, we take the neighbourhood of a barycenter of a clique as the barycenters of the other cliques it is connected to. Thus, we place additional nodes of different order in the neighbourhood of each clique in addition to the same order cliques. This information would help the cliques to have more accurate matches. The neighbours and matchings of $3$-cliques, $2$-cliques and $1$-cliques are mentioned in the Table 3 and Table 4 for both the graphs $G_{1}$ and $G_{2}$ respectively. And, the matchings shown in Figure  6, 6 and 6 are based on having the same labels for each barycenter in graph $G_{1}$ and $G_{2}$.

Appendix C Experiments

C.1 Setup

We compare the performance of our proposed method with various other matching algorithms on synthetic and real world datasets. The real world datasets are categorized in Table 5. Here, $N$ is the total number of samples with $n$ landmark points in each image to be matched. We represent random graphs on images in Figure 7 for better understanding and visualization of random graphs for our experiments. Matchings of two images for real world datasets (Table 5) are shown in Figure 14.

C.2 Effect of Affine Transformation

We created a synthetic dataset from CMU House and Hotel dataset by uniformly sampling $20\%$ and $40\%$ frames from a video sequence and performing affine transformations like rotation, reflection, scaling, and shearing. We have explained the transformations we considered for this experiment which is similar to Figure (2) and Table (1) in main paper. Table(1) in main paper shows the results on the CMU House dataset. Affine transformations on Hotel frame are shown in Figure 8. Figure 9 shows the results of matching for the remaining House (fig. 9 and 9) and Hotel synthetic dataset for all the algorithms. We observe that our method produces best results in all the cases, whereas the error for other algorithms either remains stable or increases steeply with the increase in the percentage of transformed frames in the sequence.

C.3 Effect of Occlusion

We considered two datasets with grave occlusions, mentioned in Table 5. Figures 14 and 14 show the matching of two images for both the datasets, although the matching results are shown in Table (2) in the main paper. We also created a synthetic dataset by removing $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 $20\%$ and $60\%$ of frame sequences randomly. Figures 10 and 10 show the increase in error as we remove more points from images. We also note the difference in both the results. Since we remove points from more percentage of frames in 10, there is more gradual increase in the error. This experimental setup is similar to Figure (4) in our main paper. It shows that affinity based methods like FGM and RRWM perform well but our method still consistently outperforms all the algorithms.

C.4 Effect of Frame Separation

Figures 11 and 11 show the frame separation level result of CMU House and Horse Rotate frame sequences. We select a pair of frames at a time with increase in their frame separation ($x$-axis). Here, the House dataset consists of 3D rotations of House whereas Horse Rotate dataset applies rotation with more degree of rotation as the frame separation level increases. We see that most of the algorithms performs well for both the datasets even with $0\%$ error.

C.5 Effect of k-Nearest Neighbour

In Figures 12 and 12, error and computation time of matching two frames of house are shown with different probability $p$ and nearest neighbor $k$ values. We observe that as the value of $p$ and $k$ increases, the possibility of mismatching decreases which leads to correct matching. On the other hand, the computation time increases since it increases the number of edges in the underlying graph, which in turn leads to a larger number of $d$-cliques. This also causes a marked increase in the matching algorithm’s runtime. The computation time of our algorithm considers the time of the Kuhn-Munkres algorithm, which is used as a matching algorithm to match two random clique complexes, which takes $O(n^{3})$ running time.

The overall time increases as we increase the value of $p$ and $k$, since it increases the probability of an edge occurrence between two landmark points. As the number of edges increase in a random graph, the number of $d$-cliques also increase. Due to this phenomenon, the runtime of the Kuhn-Munkres algorithm also increases.

Figure 12 shows the computation time of matching two images with varying $k$-NN for different $n$ landmark points in the image. We can clearly see that the time increases with increasing $k$ and a larger number of landmark points. Here, $60$ landmark points take maximum time for the highest value of $k$. On the other hand, if we consider lower values of $k$, even $60$ landmark points take a reasonable amount of time to match, which is comparable to lower values of $n$. Thus, we set $k$ value as low as possible for matching, depending on the complexity of the dataset.

C.6 Noise Model

We analyze the performance of our method over other pairwise algorithms for two different noise models. We follow the noise model setup mentioned in  (Feizi et al., 2016). We introduce noise in one random graph $G_{1}$ and generate a noisy version $\widetilde{G}$ to be matched with $G_{2}$. $G_{1}$ is a random graph here which is created as $G_{1}(n,p)$ with $n$ nodes and $p$ probability. We describe two noise models as follows:

Noise Model I:

[TABLE]

$\widetilde{G}$ is generated using the aforementioned equation where $A$ is a binary random symmetric matrix, whose entries are drawn from a Bernoulli distribution as $A(n,q)$ with $n$ nodes and $q$ probability and $\odot$ represents the element-wise multiplication of matrices. This model flips the node-node adjacency of $G_{1}$ with probability $q$.

Noise Model II:

[TABLE]

Again, $A$ and $B$ are binary random symmetric matrices, whose entries are drawn from the Bernoulli distribution as $A(n,q)$ and $B(n,r)$ with $n$ nodes and $q$ and $r$ probabilities, respectively. This model flips node-node adjacency of $G_{1}$ with probability $q$, and in addition it also creates edges between non-connected nodes with probability $r$.

Results of noise model I and II on CMU Hotel and Horse Shear for frame separation level is shown in Figures 13, 13 and 13, 13 respectively. We observe that our method is robust to noise for both the models as compared to other algorithms since there is a very small increase or no increase in error (%) for all the cases.