# An analysis of a random algorithm for estimating all the matchings

Jinshan Zhang<EMAIL_ADDRESS>Yan Huo
<EMAIL_ADDRESS>Fengshan Bai<EMAIL_ADDRESS>Department of Mathematical Sciences,
Tsinghua University, 100084,
Beijing, PRC.

###### Abstract

Counting the number of all the matchings on a bipartite graph has been
transformed into calculating the permanent of a matrix obtained from the
extended bipartite graph by Yan Huo, and Rasmussen presents a simple approach
(RM) to approximate the permanent, which just yields a critical ratio
O($n\omega(n)$) for almost all the 0-1 matrices, provided it’s a simple
promising practical way to compute this #P-complete problem. In this paper,
the performance of this method will be shown when it’s applied to compute all
the matchings based on that transformation. The critical ratio will be proved
to be very large with a certain probability, owning an increasing factor
larger than any polynomial of $n$ even in the sense for almost all the 0-1
matrices. Hence, RM fails to work well when counting all the matchings via
computing the permanent of the matrix. In other words, we must carefully
utilize the known methods of estimating the permanent to count all the
matchings through that transformation.

matching; permanent; critical ratio; bipartite graph; determinant; Monte-Carlo
algorithm;random algorithm; RM;fpras

###### pacs:

pacs numbers here

##  I Introduction

Let $G=(V,E)$ be a bipartite graph, where $V=V_{1}\cup V_{2}$ is the set of
vertices and $E\subset V_{1}\times V_{2}$ is the set of edges. In the
following sections we suppose $\\#V_{1}=\\#V_{2}=n$ if there’s no special
illustration. A set of edges $S\subset E$ is called a matching if no two
distinct edges $e_{1},e_{2}\in S$ contain a common vertex. $S$ is called a
k-matching if $\\#S=k$. In special case, $S$ is called a perfect matching if
$k=n$. Let $S_{k}$ be the set of k-matching in $G$ and $A(G)$ be the set of
all the k-matching, $k=0,1,\ldots,n$. For the convenience of discussion, let
$\\#S_{0}=1$, then the number of all the matchings in $G$ is
$\\#A(G)=\sum\limits_{i=0}^{n}\\#S_{k}$.

The permanent of a 0-1 $A={a_{ij},1\leq i,j\leq n}$ is defined as

$Per(A)=\sum\limits_{\pi}\prod\limits_{i=1}^{n}a_{i,\pi(i)}$ (1)

where the sum is over all the permutations $\pi$ of $[n]=\\{1,\ldots,n\\}$.
It’s well known that the permanent of an adjacent matrix of bipartite graph
equals the number of its perfect matching. Let AM(G) denote the number of all
the matchings in $G$, and $A$ be adjacent matrix of $G$. Yhuo07 has proved
that

$AM(G)=\frac{1}{n!}per\left(\begin{array}[]{cc}A&I_{n\times n}\\\ 1_{n\times
n}&1_{n\times n}\end{array}\right)$ (2)

where $I_{n\times n}$ is $n\times n$ unit matrix, $1_{n\times n}$ denotes
$n\times n$ matrix with all the elements 1. This means in order to count the
number of all the matchings of a bipartite graph with $2n$ vertices we only
need to compute the permanent of a $2n\times 2n$ corresponding matrix
transformed from adjacent matrix. The computation of permanent has a long
history and was shown to be $\\#$P-complete in Val79 . Thus, in the past $20$
years or so, many random algorithms have been developed to approximate the
permanent, which can been divided at least four categoriesChien02 : elementary
recursive algorithms(the original one is Rasmussen method(RM)) Ras94 ;
reductions to determinants Go81 ; Kar93 ; Bar99 ; Bar00 ; iterative balancing
Lin00 ; and Markov chain Monte Carlo Sin89 ; JSV01 ; KRS96 . All these methods
try to find a fully-polynomial randomized approximation scheme fpras for
computing the permanent. fpras is such a scheme which, when given
$\varepsilon$ and inputs matrix $A$, outputs a estimator(usually a unbiased
estimator)$Y$ of the permanent such that

$Pr((1-\varepsilon)per(A)\leq Y\leq(1+\varepsilon)per(A))\geq\frac{3}{4}$ (3)

and runs in polynomial time in $n$ and $\varepsilon^{-1}$, here $3/4$ may be
boosted to $1-\delta$ for any desired $\delta>0$ by running the algorithm
$O(log(\delta^{-1}))$ and taking the median of the trials Che52 . Then a
straightforward application of Chebychev’s inequality shows that running the
algorithm $O(\frac{E(Y^{2})}{E^{2}(Y)}\varepsilon^{-2})$ times and taking the
mean of the results can make the probability more than $3/4$(e.g. running
$4\frac{E(Y^{2})}{E^{2}(Y)}\varepsilon^{-2}$ times). Hence, if the critical
ratio $\frac{E(Y^{2})}{E^{2}(Y)}$ is bounded by a polynomial of inputs $A$,
we’ll get an fpras for the permanent of $A$. Another modified scheme called
fpras for almost all inputs means: choose a matrix from
$\mathcal{A}(n,1/2)$($\mathcal{A}(n,1/2)$ denotes a probability space of
$n\times n$ 0-1 matrices where each entry is chosen to be 1 or 0 with the same
probability 1/2), or equivalently choose a matrix u.a.r. from $\mathcal{A}(n)$
($\mathcal{A}(n)$ represents the set of $n\times n$ 0-1 matrices), and the
following

Pr(critical ratio of $A$ is bounded by a polynomial of the input $A$
)$=1-o(1)$ as $n\longrightarrow\infty$
holds.(Note that this is a much weaker requirement than that of an fpras). If
a proposition $P$ relating to $n$ satisfies Pr(P is true)$=1-o(n)$, we say P
holds whp(whp is the abbreviation of ”with high probability”). Thus, that
there is an fpras for almost all the matrix means the critical ratio of $A$ is
bounded by a polynomial of the input $A$ whp. A exciting result, that Markov
Chain approach led to the first fpras for the permanent of any 0-1
matrix(actually of any matrix with nonnegative entry) was shown byJSV01 .
However, its high exponent of polynomial running time makes it difficult to be
a practical method to approximate the permanent. RM and reductions to
determinants seem to be two practical approaches estimating permanent due to
their simply feasibility, and both of them have been proved to be an fpras for
almost all the 0-1 matrices. besides, Chien02 promises a good prospect on
computing permanent via clifford algebra if some difficulties can be
conquered. RM also has developed to be a kind of approaches called sequential
importance sampling way, which is widely used in statistical physics, seeBS99
.

In this paper, we’ll, by RM, compute the number of all the matchings based on
the above transformation and give its performance theoretically, say, an
analysis of critical ratio in the sense ”for almost all the 0-1 matrix” of
that matrix with a special structure. In section II, A new alternative
estimator operating directly on the adjacent matrix without any transformation
will be presented and proved to be equivalent to approximation performing on
the transformed matrix by RM. In section III, a low bound of the critical
ratio for almost all the matrices will be presented, which is larger than any
polynomial of $n$ with a certain probability. Hence, RM does not perform well
in computing the number of all the matchings as in computing the number of
perfect matching. In section IV we’ll propose some analytic results w.r.t. the
expectation and variance of the number of all the matchings of a matrix
selected u.a.r from $\mathcal{G}(m,n)$($\mathcal{G}(m,n)$ denotes the set of
bipartite graph with $\\#V_{1}=\\#V_{2}=n$ as its vertices and exact $m$
edges). These results seem likely to contribute to the upper bound of critical
ratio for almost all matrices, but the calculations are more arduous and will
be left for latter paper.

##  II An equivalent estimator

All the notations have the same meanings as those in the previous section
without special illustration. Let A an $n\times n$ 0-1 matrix be an adjacent
matrix of a bipartite graph $G=(V,E)$, $(V=V_{1}\bigcup V_{2})$. Set $Y_{A}$ a
random variable. Then RM can be stated as follows:

inputs: A an $n\times n$ 0-1 matrix;
outputs: $Y_{A}$ the estimator of permanent A;
if n=0; then

$Y_{A}=1$
else

$W=\\{j:a_{1j}=1\\}$

if $W=\emptyset$ then

$Y_{A}=0$

else

Choose $J$ u.a.r. from W

$Y_{A}=|W|Y_{1J}$

$Y_{1j}$ denotes the submatrix obtained from A by removing the 1st row and the
jth column. Note this heuristic idea comes from the Laplace’s expansion. Our
following algorithm(for easy discussion, call it AMM) is also inspired by
another expansion. we first presents our algorithm for the number of all the
matchings, and then give the explanation and proof of equivalence between AMM
and RM on the transformed matrix:

inputs: A an $n\times n$ 0-1 adjacent matrix of $G$;
outputs: $Y_{A}$ the estimator of the number of all the matchings of $G$;
if n=0; then

$Y_{A}=1$
else

$W=\\{j:a_{1j}=1\\}\bigcup\\{0\\}$

Choose $J$ u.a.r. from W

$Y_{A}=|W|Y_{1J}$

$Y_{10}$ denotes a submatrix of A by removing the 1st row(of course, it’s not
necessarily a square matrix). Define a new terminology AM on the matrix. let
$B=\\{b_{ij},1\leq i\leq m,1\leq j\leq n\\}$ an $m\times n$ matrix, $m\leq n$.
let $AM(\emptyset)=1$, by induction on $m$.

$AM(B):=AM(B_{10})+\sum\limits_{j=1}^{n}b_{1,j}B_{1j}$ (4)

Then we have the following theorem.

Theorem 1. Let A be an $n\times n$ adjacent matrix of a bipartite graph G,
Then AM(A)is the number of all the matchings of G.
Proof: It’s easy to check, when $k\geq 1$, the number of k-matching of G
equals
$\sum\limits_{i_{1},\cdots,i_{k}}\sum\limits_{\pi}a_{i_{1},\pi(i_{1})}\cdots
a_{i_{k},\pi(i_{k})}$, where $i_{1}<i_{2}\cdots<i_{k}$ chosen from
$\\{1,2,\cdots,n\\}$, $\pi$ denotes the permutation
of$\\{i_{1},i_{2},\cdots,i_{k}\\}$. Thus, the number of all the matchings
is$\sum\limits_{k=1}^{n}\sum\limits_{i_{1}<\cdots<i_{k}\\\
\subseteq\\{1,\cdots,n\\}}\sum\limits_{\pi}a_{i_{1},\pi(i_{1})}\cdots
a_{i_{k},\pi(i_{k})}+1$, where 1 denotes the number of 0-matching. Note that
if the AM(A) is written in terms of sum of elements of the matrix A, then it’s
clearly to see $AM(A)=\sum\limits_{k=1}^{n}\sum\limits_{i_{1}<\cdots<i_{k}\\\
\subseteq\\{1,\cdots,n\\}}\sum\limits_{\pi}a_{i_{1},\pi(i_{1})}\cdots
a_{i_{k},\pi(i_{k})}+1$.$\Box$

Corollary1. Let $A=\\{a_{ij}1\leq i,j\leq n\\}$ be an $n\times n$ 0-1 matrix
and $Y_{A}$ is obtained by above AMM. Then $Y_{A}$ is unbiased for AM(A),
$E(Y_{A})=AM(A)$
Proof: We prove for any $m\times n$ 0-1 matrix A, $1\leq i\leq m,1\leq j\leq
n\\}$, which will be widely used in the following proves. AMM is unbiased for
AM(A). For any fixed $n$, by induction on m, k=0,$\forall 1\leq l\leq n$,
$\forall$ a $k\times l$ 0-1 matrix A, the equation $E(Y_{A})=AM(A)$ is
trivial. Now suppose $\forall k\leq m,k\leq l\leq n$, a $k\times l$ 0-1 matrix
A has $E(Y_{A})=AM(A)$. Then when $k=m$, let $|W|=q$, we have

$\begin{split}E(Y_{A})&=\sum\limits_{j\in W}E(Y_{A}|J=j)Pr(J=j)\\\
&=\sum\limits_{j\in W}E(qY_{A_{1j}}|J=j)q^{-1}\\\ &=\sum\limits_{j\in
W}E(Y_{A_{1j}})\\\ &=\sum\limits_{j\in W}AM(A_{1j})\\\ &=AM(A).\end{split}$

$\Box$

Another simple corollary can also be obtained. To estimate the number of all
the matching in G, by RM operating on $B=\left(\begin{array}[]{cc}A&I_{n\times
n}\\\ 1_{n\times n}&1_{n\times n}\end{array}\right)$ divided by $n!$ is
equivalent to operating on A by AMM, in precise words, which can be stated as
follows.

Corollary2. Let $X_{A}$ be the output of RM operating on $A$ , $Y_{B}$ be the
output of AMM operating on transformed matrix $B$ divided by $n!$. Then
$X_{A}$ and $Y_{B}$ has the same distribution.
Proof: Note that by RM after n-th step operating on $B$,
$Y_{B}=S_{n}*Y_{1_{n\times n}}/n!$, where $S_{n}$ is a number obtained from
the first n steps, and obviously $Y_{1_{n\times n}}\equiv n!$. Hence, we have
$Y_{B}=S_{n}$. The same distribution of $S_{n}$ and $X_{A}$ can be verified
step by step.$\Box$

Corollary3. $AM(A)=\frac{1}{n!}per\left(\begin{array}[]{cc}A&I_{n\times n}\\\
1_{n\times n}&1_{n\times n}\end{array}\right)$.
Proof: This is a direct deduction of corollary2. Let $X_{A}$ be the output of
RM operating on $A$ , $Y_{B}$ be the output of AMM operating on transformed
matrix $B$ divided by $n!$.

$AM(A)=E(X_{A})=E(Y_{B})=\frac{1}{n!}per\left(\begin{array}[]{cc}A&I_{n\times
n}\\\ 1_{n\times n}&1_{n\times n}\end{array}\right)$

$\Box$
So in the following section, we’ll use AMM to compute all the matchings
instead of RM since some methodologies similar to Rasmussen can be utilized.
Another small advantage by AMM is that the critical ratio is smaller than that
directly obtained from RM. The critical ratio by RM would be $(2n)!$, see
Theorem 2.2Ras94 , while the critical ratio by AMM would be $(n+1)^{n}$.

Theorem2. Let $A=\\{a_{ij},1\leq i,j\leq n\\}$ be an $n\times n$ adjacent
matrix of a bipartite graph G, and let $X_{A}$ be the output of AMM. Then
$\frac{E(X_{A})^{2}}{E(X^{2}_{A})}\leq(n+1)^{n}$. Generally, Let A be an
$m\times n$ 0-1 matrix, $m\leq n$. $X_{A}$ be the output of AMM. Then
$\frac{E(X_{A})^{2}}{E(X^{2}_{A})}\leq(n+1)^{m}$
Proof: Induction on $m$, For any fixed $n$. $k=0$,$\forall 1\leq l\leq n$,
$\forall$ a $k\times l$ 0-1 matrix $A$, the inequation is trivial. In the case
$k=m$, let $|W|=q$, we have

$\begin{split}E(X^{2}_{A})&=\sum\limits_{j\in W}E(X^{2}_{A}|J=j)Pr(J=j)\\\
&=\sum\limits_{j\in W}E(q^{2}X^{2}_{A_{1j}}|J=j)q^{-1}\\\ &=\sum\limits_{j\in
W}E(X^{2}_{A_{1j}})q\\\ &\leq\sum\limits_{j\in
W}E(X_{A_{1j}})^{2}(n+1)^{m-1}q\\\ &\leq(\sum\limits_{j\in
W}E(X_{A_{1j}}))^{2}(n+1)^{m-1}q\\\ &=E(X_{A})^{2}(n+1)^{m}\end{split}$

$\Box$

##  III A lower bound of critical ratio for almost all the matrices

Rasmussen shows that although the critical ratio of RM is factorial in n, it
does indeed provide an fpras for almost all the matrix. However, the similar
result can not be anticipated when computing all the matchings by RM. In fact
the critical ratio for almost all the matrix would be more than
$n^{\sqrt{n}/2-1}$ with a certain probability. To prove this, we need to
define some new denotations. Since there’re two probability spaces, we use the
subscript $\sigma$ denote the calculus w.r.t. the probability space the
algorithm lies in, say, coin-tosses, and subscript $\mathcal{A}$ represent the
calculus w.r.t. the space probability the random matrices lie in.
$\mathcal{A}(m,n,p)$ denotes the probability space of all $m\times n$ 0-1
random matrices where each entry is chosen to be 1 with probability $p$,
and$\mathcal{A}(m,n)$ denotes the set of all $m\times n$ 0-1 matrices .

To obtain the mean and variance of the output of AMM on average under
probability measure $Pr_{\mathcal{A}}$, we need the following lemma.

Lemma1 Let $f(m,n)$ defined as $f(m,n)=a_{n}f(m-1,n)+c_{n}f(m-1,n-1)$, where
$m\leq n$ are two nonnegative integers, $a_{n}$ and $c_{n}$ are two infinite
positive series w.r.t. $n$. And $\forall$ $0\leq l\leq n$, $f(0,l)=1$. Then
$f(m,n)=\sum\limits_{k=1}^{m}\sum\limits_{\begin{subarray}{c}s_{0}+s_{1}+\cdots+s_{k}=m-k\\\
s_{0},\cdots s_{k}\geq 0\end{subarray}}c_{n}\cdots
c_{n-k+1}a_{n}^{s_{0}}\cdots a_{n-k}^{s_{k}}+a_{n}^{m}$
Proof: By induction on m. Obviously, the case p=0 is trivial. Suppose when
$p\leq m-1$ $\forall$ $p\leq l\leq n$,
$f(p,l)=\sum\limits_{k=1}^{p}\sum\limits_{s_{0}+s_{1}+\cdots+s_{k}=p-k}c_{l}\cdots
c_{l-k+1}a_{l}^{s_{0}}\cdots a_{l-k}^{s_{k}}+a_{l}^{p}$ holds, then when p=m,
we have

$\begin{split}a_{n}f(m-1,n)&=\sum\limits_{k=1}^{m-1}\sum\limits_{s_{0}+s_{1}+\cdots+s_{k}=m-1-k}c_{n}\cdots
c_{n-k+1}a_{n}^{s_{0}+1}\cdots a_{n-k}^{s_{k}}+a_{n}^{m}\\\
&=\sum\limits_{k=1}^{m-1}\sum\limits_{\begin{subarray}{c}s_{0}+s_{1}+\cdots+s_{k}=m-k\\\
s_{0}\geq 1\end{subarray}}c_{n}\cdots c_{n-k+1}a_{n}^{s_{0}}\cdots
a_{n-k}^{s_{k}}+a_{n}^{m}\end{split}$

and

$\begin{split}c_{n}f(m-1,n-1)&=\sum\limits_{k=1}^{m-1}\sum\limits_{s_{0}+s_{1}+\cdots+s_{k}=m-1-k}c_{n}\cdots
c_{n-k}a_{n-1}^{s_{0}}\cdots a_{n-1-k}^{s_{k}}+c_{n}a_{n-1}^{m-1}\\\
&=\sum\limits_{k=1}^{m-1}\sum\limits_{s_{1}+s_{2}+\cdots+s_{k+1}=m-1-k}c_{n}\cdots
c_{n-k}a_{n-1}^{s_{1}}\cdots a_{n-1-k}^{s_{k+1}}+c_{n}a_{n-1}^{m-1}\\\
&=\sum\limits_{k=2}^{m}\sum\limits_{s_{1}+s_{2}+\cdots+s_{k}=m-k}c_{n}\cdots
c_{n-k+1}a_{n-1}^{s_{1}}\cdots a_{n-k}^{s_{k}}+c_{n}a_{n-1}^{m-1}\\\
&=\sum\limits_{k=1}^{m}\sum\limits_{s_{1}+s_{2}+\cdots+s_{k}=m-k}c_{n}\cdots
c_{n-k+1}a_{n-1}^{s_{1}}\cdots a_{n-k}^{s_{k}}\\\
&=\sum\limits_{k=1}^{m}\sum\limits_{\begin{subarray}{c}s_{0}+s_{1}+\cdots+s_{k}=m-k\\\
s_{0}=0\end{subarray}}c_{n}\cdots c_{n-k+1}a_{n}^{s_{0}}a_{n-1}^{s_{1}}\cdots
a_{n-k}^{s_{k}}\\\ \end{split}$

From the above two equation, there holds

$\begin{split}f(m,n)&=a_{n}f(m-1,n)+c_{n}f(m-1,n-1)\\\
&=\sum\limits_{k=1}^{m}\sum\limits_{\begin{subarray}{c}s_{0}+s_{1}+\cdots+s_{k}=m-k\\\
s_{0},\cdots s_{k}\geq 0\end{subarray}}c_{n}\cdots
c_{n-k+1}a_{n}^{s_{0}}\cdots a_{n-k}^{s_{k}}+a_{n}^{m}\end{split}$

The previous $n$ can be replaced by any $l$, where $m\leq l\leq n$ $\Box$

Using lemma1 we can easily obtain two following Theorems.

Theorem3. Choose $A_{m,n}$ u.a.r. from $\mathcal{A}(m,n)$, $m\leq n$, or
equivalently let $A_{m,n}$ from $\mathcal{A}(m,n,1/2)$. Then

$E_{\mathcal{A}}(AM(A_{m,n}))=\sum\limits_{k=0}^{m}C_{m}^{k}\frac{P_{n}^{k}}{2^{k}}$

where $C_{m}^{k}=\frac{m!}{k!(m-k)!}$ and $P_{n}^{k}=\frac{n!}{(n-k)!}$
Proof: Induction on m. The case p=0, $E_{\mathcal{A}}(AM(A))=1$ is trivial.
Suppose $\forall$ $p\leq m-1$, $p\leq l\leq n$
$E_{\mathcal{A}}(AM(A_{p,l}))=\sum\limits_{k=0}^{p}C_{p}^{k}\frac{P_{l}^{p-k}}{2^{p-k}}=\sum\limits_{k=0}^{p}C_{p}^{k}\frac{P_{l}^{k}}{2^{k}}$
when p=m, $\forall$ $m\leq l\leq n$, we have

$\begin{split}E_{\mathcal{A}}(AM(A_{m,l}))&=E_{\mathcal{A}}(AM(A_{m,l}^{1,0})+\sum\limits_{j=1}^{n}a_{1,j}AM(A_{m,l}^{1,j}))\\\
&=E_{\mathcal{A}}(AM(A_{m-1,l}))+\sum\limits_{j=1}^{n}E_{\mathcal{A}}(a_{1,j})E_{\mathcal{A}}(AM(A_{m-1,l-1})\\\
&=E_{\mathcal{A}}(AM(A_{m-1,l}))+\frac{n}{2}E_{\mathcal{A}}(AM(A_{m-1,l-1})\\\
\end{split}$

Using lemma1, here $a_{l}\equiv 1$, and$c_{l}=\frac{l}{2}$ then

$\begin{split}E_{\mathcal{A}}(AM(A_{m,l}))&=\sum\limits_{k=1}^{m}\sum\limits_{\begin{subarray}{c}s_{0}+s_{1}+\cdots+s_{k}=m-k\\\
s_{0},\cdots s_{k}\geq 0\end{subarray}}c_{l}\cdots c_{l-k+1}+1\\\
&=\sum\limits_{k=1}^{m}\frac{P_{l}^{k}}{2^{k}}\sum\limits_{\begin{subarray}{c}s_{0}+s_{1}+\cdots+s_{k}=m-k\\\
s_{0},\cdots s_{k}\geq 0\end{subarray}}1+1\\\
&=\sum\limits_{k=1}^{m}\frac{P_{l}^{k}}{2^{k}}C_{m}^{k}+1\\\
&=\sum\limits_{k=0}^{m}\frac{P_{l}^{k}}{2^{k}}C_{m}^{k}\end{split}$

$\Box$
Theorem4 Choose $A_{m,n}$ u.a.r. from $\mathcal{A}(m,n)$, $m\leq n$, and let
$X_{A_{m,n}}$ be the output by AMM. Then

$E_{\mathcal{A}}(E_{\sigma}(X_{A_{m,n}}))=\sum\limits_{k=0}^{m}C_{m}^{k}\frac{P_{n}^{k}}{2^{k}}$

and

$E_{\mathcal{A}}(E_{\sigma}(X_{A_{m,n}}^{2}))=\sum\limits_{k=0}^{m}\frac{P_{n}^{k}P_{n+3}^{k}}{2^{m+k}}\sum\limits_{\begin{subarray}{c}s_{0}+s_{1}+\cdots+s_{k}=m-k\\\
s_{0},\cdots s_{k}\geq
0\end{subarray}}(n+2)^{s_{0}}(n+2-1)^{s_{1}}\cdots(n+2-k)^{s_{k}}$

Proof: The first equation is is trivial since
$E_{\sigma}(X_{A_{m,n}}^{2})=AM(A_{m,l})$. For the second one, we use
induction on m. The case p=0 is obvious. Suppose $\forall A_{p,l}$ where
$0\leq p\leq m-1$, $p\leq l\leq n$ the second equation holds. When $p=m$,
noting the fact $M=|W|-1$ is a binomial variable with parameter $l$ and
$1/2$(recall $W/\\{0\\}$ is the set of column indices with a 1 in the first
row), then

$\begin{split}E_{\mathcal{A}}(E_{\sigma}(X_{A_{m,l}}^{2}))&=\sum\limits_{q=0}^{l}E_{\mathcal{A}}(E_{\sigma}(X_{A_{m,l}}^{2})|M=q)Pr_{\mathcal{A}}(M=q)\\\
&=\sum\limits_{q=0}^{l}E_{\mathcal{A}}((q+1)\sum\limits_{j\in
W}E_{\sigma}(X_{A_{m,l}^{1j}}^{2})|M=q)Pr_{\mathcal{A}}(M=q)\\\
&=\sum\limits_{q=0}^{l}E_{\mathcal{A}}((q+1)E_{\sigma}(X_{A_{m-1,l}}^{2})+q(q+1)E_{\sigma}(X_{A_{m-1,l-1}}^{2}))Pr_{\mathcal{A}}(M=q)\\\
&=(E_{\mathcal{A}}(M)+1)E_{\mathcal{A}}(E_{\sigma}(X_{A_{m-1,l}}^{2}))+(E_{\mathcal{A}}(M^{2})+E_{\mathcal{A}}(M))E_{\mathcal{A}}(E_{\sigma}(X_{A_{m-1,l-1}}^{2}))\\\
&=(\frac{l+2}{2})E_{\mathcal{A}}(E_{\sigma}(X_{A_{m-1,l}}^{2}))+(\frac{l^{2}+3l}{4})E_{\mathcal{A}}(E_{\sigma}(X_{A_{m-1,l-1}}^{2}))\\\
\end{split}$

Using lemma1, here $a_{l}=\frac{l+2}{2}$, and $c_{l}=\frac{l^{2}+3l}{4}$.Then

$\begin{split}E_{\mathcal{A}}(E_{\sigma}(X_{A_{m,l}}^{2}))&=\sum\limits_{k=1}^{m}\frac{P_{l}^{k}P_{l+3}^{k}}{4^{k}}\sum\limits_{\begin{subarray}{c}s_{0}+s_{1}+\cdots+s_{k}=m-k\\\
s_{0},\cdots s_{k}\geq
0\end{subarray}}(\frac{l+2}{2})^{s_{0}}(\frac{l+2-1}{2})^{s_{1}}\cdots(\frac{l+2-k}{2})^{s_{k}}+(\frac{l+2}{2})^{m}\\\
&=\sum\limits_{k=0}^{m}\frac{P_{l}^{k}P_{l+3}^{k}}{2^{k+m}}\sum\limits_{\begin{subarray}{c}s_{0}+s_{1}+\cdots+s_{k}=m-k\\\
s_{0},\cdots s_{k}\geq
0\end{subarray}}(l+2)^{s_{0}}(l+2-1)^{s_{1}}\cdots(l+2-k)^{s_{k}}\\\
\end{split}$

$\Box$
Theorem5 Choose $A_{n,n}$ u.a.r. from $\mathcal{A}(n,n)$, and let
$X_{A_{n,n}}$ be the output by AMM. Then whp $h(n)\leq
E_{\mathcal{A}}(E_{\sigma}(X_{A_{n,n}}))\leq nh(n)$, where
$h(n)=\frac{(n!)^{2}}{2^{n}}\frac{2^{k^{\ast}}}{(n-k^{\ast})!(k^{\ast})^{2}}$,
$k^{\ast}=\lfloor-1+\sqrt{2n+3}\rfloor$. where $\lfloor\ast\rfloor$ denotes
the largest integer no more than $\ast$.
Proof:

$\begin{split}E_{\mathcal{A}}(E_{\sigma}(X_{A_{n,n}}))&=\sum\limits_{k=0}^{n}C_{n}^{k}\frac{P_{n}^{k}}{2^{k}}\\\
&=\sum\limits_{k=0}^{n}C_{n}^{k}\frac{P_{n}^{n-k}}{2^{n-k}}\\\
&=\frac{(n!)^{2}}{2^{n}}\sum\limits_{k=0}^{n}\frac{2^{k}}{(n-k)!(k!)^{2}}\end{split}$

and let $b_{k}=\frac{2^{k}}{(n-k)!(k!)^{2}}$, then
$\frac{b_{k}}{b_{k-1}}=\frac{2(n-k+1)}{k^{2}}$, set $\frac{b_{k}}{b_{k-1}}\geq
1$ we have $k\leq-1+\sqrt{2n+3}$, thus,
$b_{k^{\ast}}=\max\limits_{k=0,\cdots,n}b_{k}$. Thus, obviously
$\frac{(n!)^{2}}{2^{n}}b_{k^{\ast}}\leq
E_{\mathcal{A}}(E_{\sigma}(X_{A_{n,n}}))\leq
n\frac{(n!)^{2}}{2^{n}}b_{k^{\ast}}$
$\Box$
Theorem6 Choose $A_{n,n}$ u.a.r. from $\mathcal{A}(n,n)$, and let
$X_{A_{n,n}}$ be the output by AMM. Then whp

$\frac{E_{\mathcal{A}}(E_{\sigma}(X_{A_{n,n}}^{2}))}{E_{\mathcal{A}}^{2}(E_{\sigma}(X_{A_{n,n}}))}\geq
n^{(\sqrt{n}/2)}$

Proof: Numerical experiment shows the above result. however the theoretical
analysis seems so hard than until now I haven’t thought out the way to show
the comparably tight for $E_{\mathcal{A}}(E_{\sigma}(X_{A_{n,n}}^{2}))$ since
the order of
$\sum\limits_{s_{0}+s_{1}+\cdots+s_{k}=n-k}(n+2)^{s_{0}}(n+2-1)^{s_{1}}\cdots(n+2-k)^{s_{k}}$
is too difficult to gain a good lower bound. The following bound is easy to
check and the best one among methods I thought out,

$E_{\mathcal{A}}(E_{\sigma}(X_{A_{n,n}}^{2}))\geq\sum\limits_{k=0}^{n}\frac{(n!)^{2}(n+3)!}{2^{2n}}\frac{2^{k}(k+2)^{k}}{(k!)^{2}(k+3)!(n-k)!}$

However it still can’t reach the goal. Therefore, the proof of this theorem
will be left for the future.

Even if Theorem6 has been proved, unfortunately, the critical ratio for almost
all the matrices can not obtained from this theorem since two random variables
are not independent. In order to accomplish the ultimate result, we need to
calculate the $E_{\mathcal{A}}(E_{\sigma}^{2}(X_{A_{n,n}}^{2}))$. Using the
induction similar to theorem4, we can obtain the recursion of
$E_{\mathcal{A}}(E_{\sigma}^{2}(X_{A_{m,n}}^{2}))$(recall M is a binomial
variable with parameter n and $\frac{1}{2}$).

$E_{\mathcal{A}}(E_{\sigma}^{2}(X_{A_{m,n}}^{2}))=2(E_{\mathcal{A}}(M^{3})+2E_{\mathcal{A}}(M^{2})+E_{\mathcal{A}}(M))E_{\mathcal{A}}(E_{\sigma}(X_{A_{m-1,n}}^{2})E_{\sigma}(X_{A_{m-1,n-1}}^{2}))\\\
+(E_{\mathcal{A}}(M^{2})+2E_{\mathcal{A}}(M)+1)E_{\mathcal{A}}(E_{\sigma}^{2}(X_{A_{m-1,n}}^{2}))+(E_{\mathcal{A}}(M^{4})+2E_{\mathcal{A}}(M^{3})+E_{\mathcal{A}}(M^{2}))E_{\mathcal{A}}(E_{\sigma}^{2}(X_{A_{m-1,n-1}}^{2}))$

Comparing $E_{\mathcal{A}}(E_{\sigma}^{2}(X_{A_{m,n}}^{2}))$ with
$E_{\mathcal{A}}^{2}(E_{\sigma}(X_{A_{n,n}}^{2}))$ and computing their ratio
have to be done. Our main aim of doing this is to find the matrices satisfying
$E_{\sigma}(X_{A_{m,n}}^{2})\leq
E_{\mathcal{A}}(E_{\sigma}^{2}(X_{A_{m,n}}^{2}))g(n)$, where $g(n)$ is a
polynomial of $n$. However, the ratio of
$\frac{E_{\mathcal{A}}(E_{\sigma}^{2}(X_{A_{m,n}}^{2}))}{E_{\mathcal{A}}^{2}(E_{\sigma}(X_{A_{n,n}}^{2}))}$
is so large that it can’t accomplish our goal. Thus we deduce our requirement
whp to with a certain probability $p>0$, and in our results
$p=\frac{1}{2}-\varepsilon$ where $\varepsilon$ is no more than 0.02. To prove
the theorem, we need the following lemma, which will be proved in section IV.
Lemma2 Let $\mathcal{B}(m,n)$ denote the set of all $n\times n$ 0-1 matrices
with exact m 1’s, $m\gg n$. Choose $B$ u.a.r. from $\mathcal{B}(m,n)$. Then

$E(AM(B))=\sum\limits_{k=0}^{n}(C_{n}^{k})^{2}k!\frac{C_{n^{2}-k}^{m-k}}{C_{n^{2}}^{m}}$

and

$\frac{E(AM^{2}(B))}{E^{2}(AM(B))}=1+o(1),n\rightarrow\infty$

Theorem7 Choose $A_{n,n}$ u.a.r. from $\mathcal{A}(n,n)$, and let
$X_{A_{n,n}}$ be the output by AMM. Then

$Pr(\frac{E_{\sigma}(X_{A_{n,n}}^{2})}{E_{\sigma}^{2}(X_{A_{n,n}})}\geq
n^{\sqrt{n}/2-1})\geq\frac{\sum\limits_{i=(1/2+\varepsilon)n^{2}}^{n^{2}}C_{n^{2}}^{k}}{2^{n^{2}}}$

where c is a constant no more 10, and$\varepsilon\leq 0.02$.
Proof: From lemma2 we know if we set $m=(1/2+\varepsilon)n^{2}$ and
$q=\frac{C_{n^{2}-k}^{m-k}}{C_{n^{2}}^{m}}$. When n goes to infinity, noting
$k\leq n\ll m,n^{2}$, there holds

$q=\frac{C_{n^{2}-k}^{m-k}}{C_{n^{2}}^{m}}=\frac{m(m-1)\cdots(m-k)}{n^{2}(n^{2}-1)\cdots(n^{2}-k)}$

and

$\begin{split}ln(q)&=\sum\limits_{i=0}^{k-1}[ln(m-i)-ln(n^{2}-i)]\\\
&=kln(\frac{m}{n^{2}})+\sum\limits_{i=0}^{k-1}[ln(1-\frac{i}{m})-ln(1-\frac{i}{n^{2}})]\\\
&=kln(\frac{m}{n^{2}})-\sum\limits_{i=0}^{k-1}[\frac{i}{m}-\frac{i}{n^{2}}+O(\frac{i^{2}}{m^{2}})]\\\
&=kln(\frac{m}{n^{2}})-\frac{k(k-1)}{2}(\frac{1}{m}-\frac{1}{n^{2}})+O(\frac{k^{3}}{m^{2}})\end{split}$

Thus, noting that $km^{-1}\leq 2nm^{-1}=O(n^{3}m^{-2})$

$\begin{split}q&=(\frac{m}{n^{2}})^{k}exp[-\frac{k^{2}}{2}(\frac{1}{m}-\frac{1}{n^{2}})+O(\frac{n^{3}}{m^{2}})]\\\
&=(\frac{(1/2+\varepsilon)n^{2}}{n^{2}})^{k}exp[-\frac{k^{2}}{2}(\frac{1}{(1/2+\varepsilon)n^{2}}-\frac{1}{n^{2}})+O(\frac{n^{3}}{((1/2+\varepsilon)n^{2})^{2}})]\\\
&\leq e^{-1}(1/2+\varepsilon)^{k}\end{split}$

Let B selected u.a.r. from $\mathcal{B}(m,n)$ Since
$\frac{E(AM^{2}(B))}{E^{2}(AM(B))}=1+o(1)$, as $n\rightarrow\infty$
then $Pr(AM(B)<\frac{5}{6}E(AM(B)))\rightarrow 0$, as $n\rightarrow\infty$.
So, if $m\geq(1/2+\varepsilon)n^{2}$ and $\varepsilon\leq 0.02$, we have whp

$\begin{split}E_{\sigma}(X_{B}^{2})&\geq E_{\sigma}^{2}(X_{B})\\\
&=AM^{2}(B)\\\ &\geq(\frac{5}{6}E(AM(B)))^{2}\\\
&=(\frac{5}{6}\sum\limits_{k=0}^{n}(C_{n}^{k})^{2}k!\frac{C_{n^{2}-k}^{m-k}}{C_{n^{2}}^{m}})^{2}\\\
&\geq(\sum\limits_{k=0}^{n}(C_{n}^{k})^{2}k!\frac{5e^{-1}}{6}(1/2+\varepsilon)^{k})^{2}\\\
&\geq\sum\limits_{k=0}^{n}\frac{P_{n}^{k}P_{n+3}^{k}}{2^{n+k}}\sum\limits_{\begin{subarray}{c}s_{0}+s_{1}+\cdots+s_{k}=n-k\\\
s_{0},\cdots s_{k}\geq
0\end{subarray}}(n+2)^{s_{0}}(n+2-1)^{s_{1}}\cdots(n+2-k)^{s_{k}}\\\
&=E_{\mathcal{A}}(E_{\sigma}(X_{A_{n,n}})).\end{split}$

Noting
$Pr(A\in\bigcup\limits_{m\geq(1/2+\varepsilon)n^{2}}\mathcal{B}(m,n))=\frac{\sum\limits_{i=(1/2+\varepsilon)n^{2}}^{n^{2}}C_{n^{2}}^{k}}{2^{n^{2}}}$,
thus $Pr(E_{\sigma}(X_{A_{n,n}}^{2})\geq
E_{\mathcal{A}}(E_{\sigma}(X_{A_{n,n}})))\geq\frac{\sum\limits_{i=(1/2+\varepsilon)n^{2}}^{n^{2}}C_{n^{2}}^{k}}{2^{n^{2}}}$.
Using Markov’s inequality,

$Pr(E_{\sigma}(X_{A_{n,n}}^{2})\geq
nE_{\mathcal{A}}(E_{\sigma}(X_{A_{n,n}})))\leq\frac{1}{n}\rightarrow 0$

then whp $E_{\sigma}(X_{A_{n,n}}^{2})\leq
nE_{\mathcal{A}}(E_{\sigma}(X_{A_{n,n}}))$. Finally, we have

$Pr(\frac{E_{\sigma}(X_{A_{n,n}}^{2})}{E_{\sigma}^{2}(X_{A_{n,n}})}\geq\frac{1}{n}\frac{E_{\mathcal{A}}(E_{\sigma}(X_{A_{n,n}}^{2}))}{E_{\mathcal{A}}(E_{\sigma}(X_{A_{n,n}}))})\geq\frac{\sum\limits_{i=(1/2+\varepsilon)n^{2}}^{n^{2}}C_{n^{2}}^{k}}{2^{n^{2}}}$

Apply theorem6 to the above formula, we have

$Pr(\frac{E_{\sigma}(X_{A_{n,n}}^{2})}{E_{\sigma}^{2}(X_{A_{n,n}})}\geq
n^{\sqrt{n}/2-1})\geq\frac{\sum\limits_{i=(1/2+\varepsilon)n^{2}}^{n^{2}}C_{n^{2}}^{k}}{2^{n^{2}}}$

##  IV The number of all the matchings on random graph.

In this section, we consider the expectation and variance of the number of all
the matchings on G selected u.a.r. from $\mathcal{G}(m,n)$. We have the
following theorem.

Theorem8 Choose G u.a.r. from $\mathcal{G}(m,n)$, where $\mathcal{G}(m,n)$
denotes the set of bipartite graph with $\\#V_{1}=\\#V_{2}=n$ as its vertices
and exact $m$ edges, $m\gg n$, and let AM(G) denotes the number of all the
matchings in G. Then we have

$E(AM(G))=\sum\limits_{k=0}^{n}(C_{n}^{k})^{2}k!E(X_{M(k)})$

and

$E(AM^{2}(G))=\sum\limits_{k=0}^{n}\sum\limits_{i=0}^{k}(C_{n}^{k})^{2}k!\sum\limits_{p=0}^{min(i,n-k)}C_{n-k}^{p}C_{k}^{i-p}P_{n-i+p}^{p}\sum\limits_{j=0}^{i-p}C_{i-p}^{j}[F_{n-j}(i-p-j)]E(X_{M(k+i-j)})\\\
+\sum\limits_{k=1}^{n}\sum\limits_{i=0}^{k-1}(C_{n}^{k})^{2}k!\sum\limits_{p=0}^{min(i,n-k)}C_{n-k}^{p}C_{k}^{i-p}P_{n-i+p}^{p}\sum\limits_{j=0}^{i-p}C_{i-p}^{j}[F_{n-j}(i-p-j)]E(X_{M(k+i-j)})$

where $E(X_{M(k)})=C_{n^{2}-k}^{m-k}/C_{n^{2}}^{m}$ and
$F_{n}(p)=\sum\limits_{r=0}^{p}(-1)^{r}C_{p}^{r}P_{n-r}^{p-r}$
Proof: we’ll use the methodology in Jerr95 ; Let $M(k)$ be a k-matching on
$V_{1}+V_{2}$, For $G\in\mathcal{G}(m,n)$, define the random variable
$X_{M}(G)$ to be 1 if $M(k)$ is contained in G, and otherwise 0. The
expectation and second moment of $AM(G)$ is as follows.

$E(AM(G))=E(\sum\limits_{k=0}^{n}\sum\limits_{M(k)}X_{M(k)})=\sum\limits_{k=0}^{n}\sum\limits_{M(k)}E(X_{M(k)})$

and

$E(AM^{2}(G))=E((\sum\limits_{k=0}^{n}\sum\limits_{M(k)}X_{M(k)})^{2})=\sum\limits_{k=0}^{n}\sum\limits_{i=0}^{n}\sum\limits_{M(k),M^{{}^{\prime}}(i)}E(X_{M(k)}X_{M(i)}^{{}^{\prime}})$

where $\forall 0\leq k\leq n$, $M(k)$ and $M^{{}^{\prime}}(k)$ range over all
$(C_{n}^{k})^{2}k!$ k-matching’s on $V_{1}+V_{2}$. Note that

$E(X_{M(k)})=\frac{C_{n^{2}-k}^{m-k}}{C_{n^{2}}^{m}}$

The first equation follows quickly. For the second, in order to compute
$E(X_{M(k)}X_{M(i)}^{{}^{\prime}})$, we have to calculate the number of pairs
of $M(k)$ and $M^{{}^{\prime}}(i)$ as a function of the overlap
$j=|M(k)\bigcap M^{{}^{\prime}}(i)|$. For any fixed $k$, suppose $i\leq k$, we
need to compute the number of the pairs of $M(k)$ and $M^{{}^{\prime}}(i)$,
where $i=0,\cdots,k$, and $M^{{}^{\prime}}(i)$ ranges over all
$(C_{n}^{i})^{2}i!$ $i$-matching’s on $V_{1}+V_{2}$. The problem can be
equivalently stated as follows: There’re $n$ different letters and $n$
different envelopes. Among these letters, there’re exact $k(0\leq k\leq n)$
labeled letters, each of which has only one ’mother envelope’ among envelopes.
Different labeled letters have different mother envelopes. We call a $j$-fit
if there’re exact $j$ labeled letters put into its own mother envelope. Now
choose $i$($0\leq i\leq k$)letters from these $n$ letters, then put them into
$i$ envelopes, and each letter can only be put into one envelope. $\forall$
possible $j$, how many circumstances of $j$-fit are there? We can solve this
problem like this: Suppose there’re $p$ letters unlabeled and $i-p$ labeled
letters among the selected letters, obviously, $0\leq p\leq min(n-k,i)$, the
number of ways of choosing letters is $C_{n-k}^{p}C_{k}^{i-p}$. If the labeled
letters has been laid, then the number of the ways of putting $p$ unlabeled
letters is $P_{n-(i-p)}^{p}$. For any $j$($0\leq j\leq i-p$), there’re
$C_{i-p}^{j}$ ways putting exact $j$ labeled letters in its own mother
envelope. The last one we need to deal with is how many ways to put $i-p-j$
labeled letters into $n-j$ envelopes which contain all these $i-p-j$ letters’
mother envelopes, satisfying $0$-fit. By the principle of inclusion-exclusion
seeHal67 , we can easily obtain the number of the ways is $F_{n-j}(i-p-j)$,
where $F_{n}(p)=\sum\limits_{r=0}^{p}(-1)^{r}C_{p}^{r}P_{n-r}^{p-r}$. Noting
that $p$ ranges over $0$ to $min(i,n-k)$, and $j$ ranges over $0$ to $i-p$,
for each $k$ and $i\leq k$. Then

$\begin{split}\sum\limits_{M^{{}^{\prime}}(i)}E(X_{M(k)}X_{M(i)}^{{}^{\prime}})&=\sum\limits_{p=0}^{min(i,n-k)}C_{n-k}^{p}C_{k}^{i-p}P_{n-i+p}^{p}\sum\limits_{j=0}^{i-p}C_{i-p}^{j}[F_{n-j}(i-p-j)]E(X_{M(k+i-j)})\end{split}$

where $E(X_{M(k)})=C_{n^{2}-k}^{m-k}/C_{n^{2}}^{m}$ and
$F_{n}(p)=\sum\limits_{r=0}^{p}(-1)^{r}C_{p}^{r}P_{n-r}^{p-r}$.
Consider,

$\begin{split}\sum\limits_{k=0}^{n}\sum\limits_{i=0}^{n}\sum\limits_{M(k),M^{{}^{\prime}}(i)}E(X_{M(k)}X_{M(i)}^{{}^{\prime}})&=(\sum\limits_{k=0}^{n}\sum\limits_{i=0}^{k}+\sum\limits_{k=0}^{n-1}\sum\limits_{i=k+1}^{n})\sum\limits_{M(k),M^{{}^{\prime}}(i)}E(X_{M(k)}X_{M(i)}^{{}^{\prime}})\\\
&=(\sum\limits_{k=0}^{n}\sum\limits_{i=0}^{k}+\sum\limits_{k=0}^{n-1}\sum\limits_{i=k+1}^{n})\sum\limits_{M(k),M^{{}^{\prime}}(i)}E(X_{M(k)}X_{M(i)}^{{}^{\prime}})\\\
&=(\sum\limits_{k=0}^{n}\sum\limits_{i=0}^{k}+\sum\limits_{i=1}^{n}\sum\limits_{k=0}^{i-1})\sum\limits_{M(k),M^{{}^{\prime}}(i)}E(X_{M(k)}X_{M(i)}^{{}^{\prime}})\\\
&=(\sum\limits_{k=0}^{n}\sum\limits_{i=0}^{k}+\sum\limits_{k=1}^{n}\sum\limits_{i=0}^{k-1})\sum\limits_{M(k),M^{{}^{\prime}}(i)}E(X_{M(k)}X_{M(i)}^{{}^{\prime}})\\\
&=(\sum\limits_{k=0}^{n}\sum\limits_{i=0}^{k}+\sum\limits_{k=1}^{n}\sum\limits_{i=0}^{k-1})\sum\limits_{M(k)}\sum\limits_{M^{{}^{\prime}}(i)}E(X_{M(k)}X_{M(i)}^{{}^{\prime}})\\\
&=(\sum\limits_{k=0}^{n}(C_{n}^{k})^{2}k!\sum\limits_{i=0}^{k}+\sum\limits_{k=1}^{n}(C_{n}^{k})^{2}k!\sum\limits_{i=0}^{k-1})\sum\limits_{M^{{}^{\prime}}(i)}E(X_{M(k)}X_{M(i)}^{{}^{\prime}})\end{split}$

Replace $\sum\limits_{M^{{}^{\prime}}(i)}E(X_{M(k)}X_{M(i)}^{{}^{\prime}})$ by
$\sum\limits_{p=0}^{min(i,n-k)}C_{n-k}^{p}C_{k}^{i-p}P_{n-i+p}^{p}\sum\limits_{j=0}^{i-p}C_{i-p}^{j}[F_{n-j}(i-p-j)]E(X_{M(k+i-j)})$,
then the second equation is achieved.
$\Box$

Remark: To complete the proof of theorem7, we also need to know whether the
ratio $\frac{E(AM^{2}(G))}{E^{2}(AM(G))}$ goes to 1 as n goes to infinity,
adding the condition such as $m^{2}n^{-3}$ $\rightarrow$ $\infty$ as $n$
$\rightarrow$ $\infty$. We guess such a result is right, however the calculus
seems very difficult. And this result also contributes to the upper bound of
critical ratio for almost all the matrices.

###### Acknowledgements.

## References

  * (1) M. Hall JR.  Combinatorial Theory,  Blaisdell, Waltham Massachusetts (1967).
  * (2) Valiant.  The complexity of computing the permanent,  _Theoretical Computer Science_. 8, 189-201 (1979).
  * (3) S. Chien, L. Rasmussen and A. Sinclair.  Clifford Algebras and approximating the permanent,  _Proceedings of the 34th Annual Symposium on Theory of Computing (STOC)_. 222 C231 (2002).
  * (4) L. E. Rasmussen.  Approximating the Permanent: a Simple Approach,  _Random Structures and Algorithms_. 5, 349-361 (1994).
  * (5) C. Godsil and I. Gutman.  On the matching polynomial of a graph,  _Algebraic Methods in Graph Theory_. 241-249 (1981).
  * (6) A. Frieze and M. Jerrum.  An Analysis of A Monte Carlo Algorithm For Estimating the Permanent,  _Combinatorica_. 15(1), 67-83 (1995).
  * (7) N. Karmarkar, R. Karp, R. Lipton, L. Lovász, and M. Luby,.  A Monte-Carlo algorithm for estimating the permanent,  _SIAM Journal on Computing_. 22, 284-293 (1993).
  * (8) Yan Huo, Heng Liang, Siqi Liu and Fengshan Bai.  Approximating the monomer-dimer constants through matrix permanent,  _arXiv:0708.1641v2:cond-mat.stat-mech_. (2007).
  * (9) A. Barvinok.  Polynomial time algorithms to approximate permanents and mixed discriminants within a simply exponential factor,  _Random Struct. Algorithms_. 14, 29-61 (1999).
  * (10) H. Chernoff.  A measure ofasymptotic efficiencyfor tests ofa hypothesis based on the sum ofobservations,  _Ann. Math. Stat._ 23, 493-509 (1952).
  * (11) A. Barvinok.  New permanent estimators via non-commutative determinants,  Preprint. (2000).
  * (12) N. Linial, A. Samorodnitsky and A. Wigderson.  A deterministic strongly polynomial algorithm for matrix scaling and approximate permanents,  _Combinatorica_. 20, 545-568 (2000).
  * (13) M. Jerrum and A. Sinclair.  Approximating the permanent, ,  _SIAM Journal on Computing_. 18, 1149-1178 (1989).
  * (14) I.Beichl and F. Sullivan.  Approximating the Permanent via Importance Sampling with Application to the Dimer Covering Problem,  _Journal of Computational Physics_. 149, 128-147 (1999).
  * (15) Shmuel Friedland and Daniel Levy.  A polynomial-time approximation algorithm for the number of k-matchings in bipartite graphs,  _arXiv:cs.CC/0607135 v1_. 392-401 (2006).
  * (16) M. Jerrum, A. Sinclair and E. Vigoda.  A polynomial-time approximation algorithm for the permanent of a matrix with non-negative entries,  _Proceedings of the 33rd ACM Symposium on Theory of Computing_. 712 C721 (2001).
  * (17) A. Frieze and S. Suen.  Counting the number of hamilton cycles in random digraphs,  _Random Structrues and Algorithms_.(1992).
  * (18) B. Bollobás.  Random Graph,  _Cambridge University Press_. Second Edition (2001).
  * (19) C. Kenyon, D. Randall, and A. Sinclair.  Approximating the number of dimer coverings of a lattice,  _Journal of Statistical Physics_. 83, 637-659 (1996).
  * (20) George E. Andrews, Richard Askey and Ranjan Roy.  Special Function,  _Cambridge University Press_. (1999).


#  CLIPPER+: A Fast Maximal Clique Algorithm
for Robust Global Registration

Kaveh Fathian1, Tyler Summers2 Manuscript received: September 29, 2023;
Revised December 28, 2023; Accepted January 15, 2024.This paper was
recommended for publication by Editor Javier Civera upon evaluation of the
Associate Editor and Reviewers’ comments. This work is supported by the United
States Air Force Office of Scientific Research under Grant FA9550-23-1-0424,
by the National Science Foundation under Grant ECCS-2047040, and by the ARL
SARA under Grant W911NF2420029.1K. Fathian is with the Department of Computer
Science, Colorado School of Mines. Email<EMAIL_ADDRESS>Summers
is with the Department of Mechanical Engineering, University of Texas at
Dallas. Email<EMAIL_ADDRESS>
###### Abstract

We present CLIPPER+, an algorithm for finding maximal cliques in unweighted
graphs for outlier-robust global registration. The registration problem can be
formulated as a graph and solved by finding its maximum clique. This
formulation leads to extreme robustness to outliers; however, finding the
maximum clique is an NP-hard problem, and therefore approximation is required
in practice for large-size problems. The performance of an approximation
algorithm is evaluated by its computational complexity (the lower the runtime,
the better) and solution accuracy (how close the solution is to the maximum
clique). Accordingly, the main contribution of CLIPPER+ is outperforming the
state-of-the-art in accuracy while maintaining a relatively low runtime.
CLIPPER+ builds on prior work (CLIPPER [1] and PMC [2]) and prunes the graph
by removing vertices that have a small core number and cannot be a part of the
maximum clique. This will result in a smaller graph, on which the maximum
clique can be estimated considerably faster. We evaluate the performance of
CLIPPER+ on standard graph benchmarks, as well as synthetic and real-world
point cloud registration problems. These evaluations demonstrate that CLIPPER+
has the highest accuracy and can register point clouds in scenarios where over
$99\%$ of associations are outliers. Our code and evaluation benchmarks are
released at https://github.com/ariarobotics/clipperp.

###### Index Terms:

Localization; SLAM; RGB-D Perception

##  I Introduction

Data association is broadly defined as the correspondence of identical/similar
elements across sets of data, and is a key component of many robotics and
computer vision applications, such as localization and mapping [3, 4, 5],
point cloud registration [6, 7], shape alignment [8], object detection [9],
data fusion [10, 11], and multi-object tracking [12]. In these applications,
it is crucial that data association is solved correctly and fast.

In point cloud registration, for example, we seek to find the rigid
transformation (rotation/translation) that aligns two sets of 3D points. This
requires associating points in one set with their corresponding points in the
other set. Local registration techniques such as the Iterative Closest Point
(ICP) algorithm [13] associate points based on their nearest neighbor. These
associations are generally wrong if the point clouds are not aligned well
initially, leading to wrong registration. Better associations can be
established by matching descriptors that are computed around each point in the
point cloud and describe the local geometry and appearance of a point (e.g.,
classical FPFH [14] or modern learning-based 3DMatch [15]). However, due to
noise, repetitive patterns, small overlap between the point clouds, etc.,
these putative associations can have extreme outlier ratios (e.g., FPFH
associations in Section V-C are $99\%$ outliers/wrong). In these high-outlier
regimes, existing outlier rejection techniques (such as the general RANSAC
framework [16] or specific frameworks for point cloud registration [17])
either return wrong results or have impractical runtime (e.g., RANSAC’s
runtime grows exponentially in outlier ratio [6]).

To address these issues, we present the CLIPPER+ algorithm. CLIPPER+
formulates the data association problem as a graph, in which the
inlier/correct associations are the maximum clique. This formulation is robust
to high outlier ratios and applicable to any problem that admits invariants
(see Section III) such as point, line, and plane could registration [1, 18,
19]. To address the high computational complexity of finding the maximum
clique (NP-hardness), CLIPPER+ finds an approximate solution instead, which is
obtained from combining an improved version of our prior work, CLIPPER [1],
and the greedy maximal clique algorithm in [2]. CLIPPER+ runs in polynomial
time and outperforms state-of-the-art algorithms in maximum clique estimation
accuracy (Section V-A). Further, CLIPPER+ solutions are shown to be exact
(i.e., the maximum clique) in over $99\%$ of the point cloud registration
trials (Section V-C).

Figure 1: An example of maximum clique formulation for robust global
registration. (Left) Putative associations (red: outliers, green: inliers).
(Right) Graph formulation with maximum clique indicating inlier associations.

Contributions: In summary, this work’s contributions are:

  * •

An improved solver (Algorithm 2) leading to higher accuracy over our prior
work CLIPPER [1].

  * •

The new CLIPPER+ algorithm as the combination of the greedy algorithm in [2]
and CLIPPER [1], further improving both runtime and accuracy.

  * •

Evaluations demonstrating superior accuracy of CLIPPER+ over state-of-the-art
on maximum clique problems, and correct registration of real-world point
clouds in regimes with over $99\%$ outliers.

  * •

Efficient C++ implementation of all algorithms (open-source code will be
released).

##  II Related Work

Robust registration: In robotic applications of data association, such as
registration, formulations that are robust to uncertainties (such as extreme
noise and outliers) are based on globally optimal solvers, such as Branch and
Bound [20, 21, 6], and combinatorial approaches for maximum consensus [22].
These methods can tolerate extreme outlier ratios, but have worst-case
exponential complexity and are slow in practice. To address this issue, fast
heuristics such as RANSAC [16] and its variants [23], graduated non-convexity
[17, 24, 25], or iterative local solvers for M-estimation [26, 27] are used,
however, they are prone to failures in high-outlier regimes [6, 7, 28].

Graph-theoretic formulation: The dichotomy between robustness and runtime can
be resolved in some data association settings that admit “invariants” (see
Section III, and [19, 1, 18]) by formulation in a graph-theoretic framework.
Recent algorithms that leverage this framework are CLIPPER [1], TEASER [7],
and ROBIN [19] for global point cloud registration, robust to more than $99\%$
outliers and suitable for real-time applications, and PCM [4] for pairwise-
consistent loop closures in SLAM. This graph-theoretic formulation can be
traced back to Bailey et al. [29] for 2D LiDAR registration as a maximum
common subgraph problem, for which the maximum clique indicated the correct
data association. Leordeanu and Hebert [30] extended this graph-theoretic
framework to weighted consistency graphs. Enqvist et al. [31] noted the
suboptimality of [30] and proposed a vertex covering formulation, essentially
an alternative to the maximum clique formulation of [29]. Parra et al. [32]
proposed a practical maximum clique algorithm for geometric consistency based
on branch and bound and graph coloring. Recently, Zhang et al. [33] have shown
that the registration problem can be solved by relaxing the maximum clique
constraint and using the maximal cliques instead.

Maximum clique estimation: Finding the maximum clique is a well-known NP-hard
problem in its full generality [34], meaning that the complexity of finding
the solution using exact algorithms (e.g., [2]) grows exponentially in graph
size, which is impractical for data association applications. An approximate
(maximal clique) solution can however be found by polynomial-time algorithms.
Notable algorithms in this setting are Pelillo [35] and Ding et al. [36],
based on the Motzkin-Straus formulation [37], Belachew and Gillis [38] based
on symmetric rank-one nonnegative Matrix approximation, and our prior CLIPPER
algorithm [1] based on a continuous relaxation.

##  III Graph-Theoretic Robust Registration

The key idea for creating robustness to extreme outlier ratios is to find the
largest set of jointly consistent data and/or associations. This problem can
be formulated as a graph, in which this set is represented by the maximum
clique. In what follows, we introduce this graph-theoretic framework and its
use case for point cloud registration.

Point cloud registration: The objective of point cloud registration is to find
the rigid transformation that aligns a set of points to their corresponding
points in another set. The main challenge is finding the correct
correspondences as usually only a subset of the points match, and, the points
do not align perfectly due to noise. An outlier can be a point in one set that
does not have a counterpart in the other set, or an association that matches
wrong points across the sets. Our focus here is associations, and henceforth
outliers imply outlier associations. Fig. 1 illustrates a point cloud
registration example where we seek to find the blue bunny in the cluttered
point cloud.

Maximum likelihood solution: In the absence of prior knowledge and when
outliers are random, unbiased, and unstructured, the largest set of jointly
consistent associations are inliers for the maximum likelihood solution. For
point cloud registration, two associations are consistent if their endpoints
are equidistant, and therefore can be aligned by a rigid transformation. In
Fig 1, the green associations align $3$ points while the red associations
align $2$ points. Thus, the green associations are the inliers and can be used
to compute the correct maximum likelihood solution.

Graph formulation: Finding the largest set of jointly consistent associations
can be formulated as a maximum clique problem. Given $n$ associations, the
consistency graph is a graph of $n$ vertices where each vertex represents an
association. An edge between two vertices indicates that the associations are
consistent. In Fig. 1, an edge between two vertices of the consistency graph
(shown on the right) indicates consistency of corresponding associations. For
example, there is an edge between vertices/associations $2$ and $3$ as their
endpoints are equidistant ($d=d^{\prime}$), while there is no edge between
vertices/associations $1$ and $3$. Given two associations with endpoint
distances $d$ and $d^{\prime}$, due to noise, often a consistency threshold
$\epsilon$ is used where if $|d-d^{\prime}|<\epsilon$, the associations are
deemed consistent.

A clique is a subset of vertices where every pair of vertices within that
subset is connected by an edge. The maximum clique is the clique with the
largest number of vertices. A maximal clique is a clique that is not contained
in a larger clique. The maximum clique in Fig. 1 consists of vertices $2$,
$3$, and $5$. Vertices $1$ and $4$ form a maximal clique.

Invariants: The graph-theoretic framework can be applied to a broad array of
data association problems in robotics. An invariant is a quantity that remains
unchanged across the sets. The invariant used for the point cloud registration
example above was the Euclidean distance between the endpoints. Invariants can
be defined to register lines [1], planes [18], 2D-3D visual features [19],
etc. Examples in this work, however, focus on point cloud registration.

##  IV Maximal Clique Algorithms

Algorithm 1 Degeneracy-Ordered Greedy Maximal Clique

1:Input $A$: adjacency matrix; $K$: core numbers

2:Output $C$: vertices that form maximal clique

3:$C\leftarrow\\{\\}$; $c_{\max}\leftarrow 0$

4:% Sort vertices by core number in descending order:

5:$(v_{1},\dots,v_{n})\leftarrow$ sort vertices $v_{i}$ s.t. $K(v_{i})\geq
K(v_{j})$ if $i<j$

6:for $i=1:n$ do

7: if $K(v_{i})\geq c_{\max}$ then

8: % Neighbors of $v_{i}$ with core numbers $\geq c_{\max}$:

9: $S\leftarrow\\{v_{j}:A(v_{i},v_{j})=1,K(v_{j})\geq c_{\max}\\}$

10: $(s_{1},\dots,s_{k})\leftarrow$ sort $s_{i}\in S$ descending by core
number

11: $C^{\prime}\leftarrow\\{\\}$

12: for $j=1:k$ do % For each vertex in sorted $S$

13: if $A(c_{i},s_{j})=1~{}\forall_{c_{i}\in C\textquoteright}$ then %
$C\textquoteright\cup\\{s_{j}\\}$ is clique

14: $C\textquoteright\leftarrow\\{C\textquoteright,s_{j}\\}$ % Add $s_{j}$ to
$C\textquoteright$

15: if $|C\textquoteright|>c_{\max}$ then % Found a larger clique

16: $C\leftarrow C\textquoteright$; $c_{\max}\leftarrow|C^{\prime}|$ % Update
the output

We present the main algorithmic contributions of this work in subsections IV-B
and IV-C, where we discuss the maximal clique algorithms based on a continuous
relaxation and CLIPPER+. Subsection IV-A is a review of the approach in [2],
and is not an algorithmic contribution. However, our C++ implementation of
this algorithm improved the performance and accuracy compared to its original
implementation.

###  IV-A Degeneracy-Ordered Greedy Maximal Clique Algorithm

Algorithm 1 presents a greedy approach for finding a maximal clique. Starting
from an empty clique (line 3), we grow the clique one vertex at a time by
looping through the vertices (line 6). For each vertex $v$ that this loop
examines, the algorithm adds $v$ to the current clique if it is connected to
every vertex that is already in the clique, and discards $v$ otherwise (lines
13-16). This algorithm has the overall runtime of $\mathrm{O}(\delta|E|)$,
where $\delta$ is the maximum vertex degree and $|E|$ is the number of edges
[2].

The maximal clique returned by the greedy approach depends on the initial
vertex chosen to grow the clique, and the ordering of vertices (as they are
sequentially examined to be added to the current clique or discarded). A
descending ordering of vertices by their core number greatly improves the odds
of finding a large maximal clique (and possibly the maximum clique), as
leveraged in Algorithm 1 (lines 5 and 10). Mathematically, the core number or
degeneracy of a vertex $v$ is the largest integer $k$ such that the degree of
$v$ remains non-zero when all vertices of degree less than $k$ are recursively
removed from the graph. The core number of vertices can be computed
efficiently in $\mathrm{O}(|E|)$ by Batagelj and Zaversnik’s algorithm [39].

###  IV-B Continuous-Relaxation Maximal Clique Algorithm

Algorithm 2 Continuous-Relaxation Maximal Clique

1:Input $A$: adjacency matrix; $\bar{u}$: initial guess

2:Output $C$: vertices that form maximal clique

3:Parameters $\sigma\leftarrow 0.01$; $\beta\leftarrow 0.5$; $tol\leftarrow
10^{-8}$

4:$M\leftarrow A+I$ % Add identity matrix $I$

5:$\bar{M}\leftarrow\mathbf{1}-M$ % Binary complement of $M$

6:$u\leftarrow\max(\bar{u}/\|\bar{u}\|,0)$; $d\leftarrow d_{0}$;
$\alpha\leftarrow 1$

7:while $u$ not in binary state do

8: $M_{d}\leftarrow M-d\,\bar{M}$

9: $F\leftarrow u^{\top}M_{d}u$

10: % Gradient projected on $S^{n}$ tangent bundle:

11: $\nabla F_{\perp}(u)=2(I-uu^{\top})M_{d}u$

12: $\Delta u\leftarrow tol$; $\Delta F\leftarrow tol$

13: while $\Delta u\nless tol$ or $\Delta F\nless tol$ do

14: $Armijo\leftarrow\mathrm{False}$

15: while $Armijo=\mathrm{False}$ do % Backtracking line search

16: $u_{+}\leftarrow u+\alpha\nabla F_{\perp}(u)$ % Solution update

17: $u_{+}\leftarrow\max(u_{+}/\|u_{+}\|,0)$ % Retract to
$\mathbb{R}^{n}_{+}\cap S^{n}$

18: $F_{+}\leftarrow u_{+}^{\top}M_{d}u_{+}$

19: $\nabla F_{\perp}(u_{+})=2(I-u_{+}u_{+}^{\top})M_{d}u_{+}$

20: $\Delta F\leftarrow F_{+}-F$; $\Delta u\leftarrow u_{+}-u$

21: $Armijo\leftarrow(\Delta F\geq\sigma\nabla F_{\perp}(u)^{\top}\Delta u)$

22: if $Armijo=\mathrm{False}$ then

23: $\alpha\leftarrow\alpha\,\beta$ % Reduce $\alpha$

24: else $\alpha\leftarrow\alpha/\sqrt{\beta}$ % Increase $\alpha$

25: $u\leftarrow u_{+}$; $\nabla F_{\perp}(u)\leftarrow\nabla
F_{\perp}(u_{+})$; $F\leftarrow F_{+}$

26: $d\leftarrow d+\Delta d$ % Increase $d$

27:$C\leftarrow\\{i:u_{i}>0\\}$ % Vertices that form maximal clique

Optimization formulation: The maximum clique problem in an undirected and
unweighted graph of $n$ vertices can be formulated as

$\displaystyle\begin{array}[]{ll}\underset{u\in\\{0,1\\}^{n}}{\text{maximize}}&\sum_{i=1}^{n}{u_{i}}\\\
\text{subject
to}&u_{i}\,u_{j}=0\quad\text{if}~{}A(i,j)=0,~{}\forall_{i,j}\end{array}$ (3)

where $A\in\\{0,1\\}^{n\times n}$ is the adjacency matrix with $A(i,j)=1$ if
and only if vertices $u_{i}$ and $u_{j}$ are connected. The optimization
variable $u$ is a binary vector of $n$ elements, where $1$ entries indicate
vertices that form a clique. As $u$ is binary, the constraint $u_{i}\,u_{j}=0$
if $A(i,j)=0$ implies that if vertices $u_{i}$ or $u_{j}$ are disconnected,
then at most one of them can be selected. For example, consider the graph in
Fig. 1 with the adjacency matrix $A$ and solution candidates $u,\bar{u}$ as

$A=\begin{bmatrix}0&0&0&1&0\\\ 0&0&1&0&1\\\ 0&1&0&0&1\\\ 1&0&0&0&0\\\
0&1&1&0&0\\\ \end{bmatrix},~{}u=\begin{bmatrix}0\\\ 1\\\ 1\\\ 0\\\ 1\\\
\end{bmatrix},~{}\bar{u}=\begin{bmatrix}1\\\ 0\\\ 0\\\ 1\\\ 0\\\
\end{bmatrix}.$ (4)

Both $u$ and $\bar{u}$ satisfy the constraints in (3). Solution $u$ is the
global optimum with the objective value of $3$ (the size of the maximum
clique, and the number of $1$ entries in $u$), while $\bar{u}$ gives the
objective value of $2$.

By defining $M\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny def}}}{=}}A+I$,
where $I$ is the identity matrix of appropriate size, it is straightforward to
show that problem (3) is equivalent to

$\displaystyle\begin{array}[]{ll}\underset{u\in\\{0,1\\}^{n}}{\text{maximize}}&\dfrac{u^{\top}M\,u}{u^{\top}u}\\\
\text{subject
to}&u_{i}\,u_{j}=0\quad\text{if}~{}M(i,j)=0,~{}\forall_{i,j}.\end{array}$ (7)

This follows by observing that in (3) and (7) the constraints are identical in
the sense that disconnected vertices cannot be selected jointly in the
solution. Further, if $u^{*}$ is an optimal solution of (7) and it has $m$ one
entries, then the objective values of (7) and (3) will both be identical and
equal to $m$ (e.g., $u$ and $\bar{u}$ in (4) give objective values of $3$ and
$2$ in (7)).

Continuous relaxation: To overcome the NP-hardness of problem (7), an
approximate solution can be obtained by a continuous relaxation, where the
binary domain ${u\in\\{0,1\\}^{n}}$ is relaxed to the set of non-negative real
numbers ${u\in\mathbb{R}^{n}_{+}}$. Gradient-based optimization routines with
polynomial time complexity can solve the relaxed problem, and the relaxed
solution can be projected/rounded to the nearest binary. The issue of this
approach is that there is no guarantee that the binarized solution is a
feasible solution that satisfies the constraints of the original problem (7).

Our proposed relaxation of (7) that addresses this issue is

$\displaystyle\begin{array}[]{ll}\underset{u\in\mathbb{R}^{n}_{+}}{\text{maximize}}&\dfrac{u^{\top}M_{d}\,u}{u^{\top}u},\end{array}$
(9) $\displaystyle M_{d}(i,j)\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}\left\\{\begin{array}[]{ll}M(i,j)&\text{if}~{}~{}M(i,j)\neq 0\\\
-d&\text{if}~{}~{}M(i,j)=0\end{array}\right.$ (12)

where $d>0$ is a positive scalar. Problem (9) can be written equivalently as

$\displaystyle\begin{array}[]{ll}\underset{u\in\mathbb{R}^{n}_{+}}{\text{maximize}}&F(u)\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}u^{\top}M_{d}\,u\\\ \text{subject to}&\|u\|=1\end{array}$ (15)

where $\|\cdot\|$ is the $\ell_{2}$ vector norm.111Vector $u$ can be written
as $u=cu^{\prime}$, where $c=\|u\|$, and $u^{\prime}=u/c$ is a unit-norm
vector. Replacing $u=cu^{\prime}$ in (9) gives
$\frac{u^{\top}M_{d}u}{u^{\top}u}=\frac{c^{2}{u^{\prime}}^{\top}M_{d}u^{\prime}}{c^{2}{u^{\prime}}^{\top}u^{\prime}}={u^{\prime}}^{\top}M_{d}\,u^{\prime}$,
demonstrating the equivalency to (15).

Our relaxation (15) is inspired by Belachew and Gillis [38], which integrates
the constraints in (7) into (15) using the matrix $M_{d}$. The difference of
this work with [38] is that in [38] the relaxation
$\min_{u\in\mathbb{R}_{+}^{n}}\|M_{d}-uu^{\top}\|_{F}^{2}$ is used.
Intuitively, any entry $M_{d}(i,j)=-d$ penalizes joint selection of
disconnected vertices $u_{i}$ and $u_{j}$ by the amount $-2\,d$ in the
objective. Hence, as $d$ increases, the entries of $u$ that violate clique
constraints converge to zero.

Optimality guarantees: When $d\geq n$ (see (12)), the optima of the proposed
relaxation (15) are theoretically guaranteed to correspond to the optima of
the original maximum clique problem (7). That is, a local optimum of (15)
corresponds to a maximal clique, and a global optimum of (15) corresponds to a
maximum clique. It is by no means trivial to prove this statement, and the
interested reader should refer to Theorems 3-5 in [38] that prove this
statement for the relaxation
$\min_{u\in\mathbb{R}_{+}^{n}}\|M_{d}-uu^{\top}\|_{F}^{2}$, and Theorem 2 in
[38] that establishes connections to our proposed relaxation (15).

It is interesting to point out that any solution of (15) has a binary state
(this follows from the analysis in [38]). That is, the entries of a solution
vector $u^{*}\in\mathbb{R}^{n}_{+}$ are either $0$ or equal to a positive
scalar $c>0$. These positive entries are indicators of vertices that form a
maximal clique. Noting that (15) can be (locally) solved in polynomial time by
using a gradient-based solver, and the one-to-one correspondence between the
optima of (15) and the (NP-hard) maximum clique problem (7), one may think
that the maximum clique problem can be solved in polynomial time. Note,
however, that (15) is a nonlinear optimization problem, potentially with many
local optima, and depending on the initial guess the solver can converge to a
local optimum instead of the global optimum. Finding the global optimum of
(15), hence, remains an NP-hard problem.

Algorithm 3 CLIPPER+

1:Input $A$: adjacency matrix

2:Output $C$: vertices that form maximal clique

3:$K\leftarrow$ compute core number of vertices % From [39]

4:$C\leftarrow$ Algorithm1($A,K$) % Degeneracy-ordered greedy clique

5:$\mathcal{I}\leftarrow\\{i:K(v_{i})\geq|C|\\}$ % Vertices with core number
$\geq$ greedy clique size

6:$A\textquoteright\leftarrow A(\mathcal{I},\mathcal{I})$ % Prune graph (only
keep vertices in $\mathcal{I}$)

7:if $A\textquoteright=\emptyset$ then Terminate % Maximum clique found

8:% Binary complement of greedy clique in pruned graph:

9:$\bar{u}\leftarrow(\bar{u}_{i}:\bar{u}_{i}=0~{}\mathrm{if}~{}i\in
C,~{}\mathrm{else},~{}\bar{u}_{i}=1,~{}\forall_{i\in\mathcal{I}})$

10:$C^{\prime}\leftarrow$ Algorithm2($A^{\prime},\bar{u}$) % Continuous-
relaxation clique

11:if |C’| > |C| then $C\leftarrow C^{\prime}$ % Return larger clique

TABLE I:  Comparisons of the maximum clique estimation accuracy and runtime
for CLIPPER+ on DIMACS benchmark [40] (bold is best).

Benchmark |  |  |  |  | Pelillo |  | Ding |  | Belachew |  | CLIPPER |  | Greedy |  | Optim |  | CLIPPER+
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Graph | $n$ | $s$ | $\omega_{\mathrm{gt}}$ |  | $t\downarrow$ | $r\uparrow$ |  | $t\downarrow$ | $r\uparrow$ |  | $t\downarrow$ | $r\uparrow$ |  | $t\downarrow$ | $r\uparrow$ |  | $t\downarrow$ | $r\uparrow$ |  | $t\downarrow$ | $r\uparrow$ |  | $t\downarrow$ | $r\uparrow$
`C125.9` | 125 | 0.1 | 34 |  | 4.1 | 0.97 |  | 4.5 | 0.56 |  | 10.9 | 1 |  | 9.1 | 1 |  | 0.7 | 0.85 |  | 0.8 | 1 |  | 1.6 | 1
`C250.9` | 250 | 0.1 | 44 |  | 4.7 | 0.89 |  | 3.6 | 0.34 |  | 7.8 | 0.95 |  | 6.8 | 0.95 |  | 2.7 | 0.8 |  | 6.4 | 0.95 |  | 8.7 | 0.95
`brock200_2` | 200 | 0.5 | 12 |  | 2.2 | 0.67 |  | 51.8 | 0.75 |  | 4.3 | 0.83 |  | 4 | 0.83 |  | 1.4 | 0.83 |  | 2.6 | 0.83 |  | 3.8 | 0.83
`brock200_4` | 200 | 0.34 | 17 |  | 1.7 | 0.76 |  | 9.4 | 0.88 |  | 4.2 | 0.88 |  | 8.9 | 0.65 |  | 1.7 | 0.82 |  | 29.4 | 0.82 |  | 3.6 | 0.94
`gen200_p0.9_44` | 200 | 0.1 | 44 |  | 6.1 | 0.75 |  | 3.9 | 0.2 |  | 10.4 | 0.84 |  | 10.9 | 0.89 |  | 2.2 | 0.73 |  | 4.7 | 0.89 |  | 6.4 | 0.89
`gen200_p0.9_55` | 200 | 0.1 | 55 |  | 4.6 | 0.69 |  | 2.7 | 0.38 |  | 4.9 | 0.69 |  | 6.8 | 1 |  | 2.0 | 0.64 |  | 1.6 | 1 |  | 3.7 | 1
`keller4` | 171 | 0.35 | 11 |  | 2.3 | 0.64 |  | 4.0 | 0.64 |  | 2.8 | 0.73 |  | 2.8 | 0.64 |  | 1.0 | 0.82 |  | 35.4 | 0.64 |  | 5.8 | 0.82
`p_hat300-1` | 300 | 0.76 | 8 |  | 1.3 | 0.75 |  | 33.9 | 0.88 |  | 6.5 | 1 |  | 4.3 | 1 |  | 2.1 | 0.88 |  | 4.9 | 1 |  | 10.7 | 1
`p_hat300-2` | 300 | 0.51 | 25 |  | 2.8 | 0.96 |  | 11.3 | 0.92 |  | 10.9 | 1 |  | 15.9 | 0.96 |  | 3.7 | 0.84 |  | 6.5 | 1 |  | 10.1 | 1

$n$: Number of graph vertices. $s$: Graph sparsity. $\omega_{\mathrm{gt}}$:
Ground truth maximum clique size. $t$: Runtime (milliseconds); the lower, the
better. $r$: Maximum clique accuracy ratio
($\hat{\omega}/\omega_{\mathrm{gt}}$); the closer to $1$, the better.

Optimization algorithm: We present a custom solver for (15) based on a
projected gradient ascent approach, as described in Algorithm 2. Our approach
is similar to the algorithm in our previous work CLIPPER [1]; however,
compared to [1], Algorithm 2 is significantly improved and uses the Armijo
procedure for selecting the appropriate step size, which leads to more
accurate results (as we will show in the comparisons).

Problem (15) is nonlinear with many local optima in general (corresponding to
maximal cliques). To improve the odds of finding the global optimum (maximum
clique) and escape local optima, in Algorithm 2 we use a homotopy approach
where we increase $d$ incrementally in an outerloop (lines 7-26). As the
penalty parameter $d$ increases incrementally by $\Delta d$ in each iteration
of the outerloop (line 26), the elements of $u$ that violate the clique
constraints are penalized further and $u$ converges to a feasible solution.
This process continues until $d$ is large enough ($d\geq n$) and $u$ converges
to a binary state, corresponding to a maximal clique. The $\Delta d$
increments can be chosen as a small constant (as done in [38]), however, in
our implementation, we use a greedy scheme where we increase $d$ until the
smallest element of $u$ that violates the clique constraint goes to zero in
the next iteration. In the final step (line 27), the vertices of the maximal
clique are identified as the non-zero elements of $u$.

The innerloop (lines 13-25) ensures that for each $d$ increment enough
iterations of the gradient ascent are performed for the solution $u$ to reach
a steady state. Noting that the constraint manifold of the optimization
problem (15) is $\mathbb{R}^{n}_{+}\cap S^{n}$, where $S^{n}$ is the unit
sphere, to speed up the convergence, we project the gradient $\nabla
F(u)\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny def}}}{=}}2M_{d}\,u$ onto
the tangent bundle of $S^{n}$ at $u$ and move along the orthogonal projection
$\nabla F_{\perp}(u)=2(I-uu^{\top})M_{d}u$ (line 11).

To find an appropriate step size $\alpha$ along the projected gradient, we use
backtracking line search (lines 15-24) with the Armijo procedure (lines
21-24), which guarantees a sufficient increase in the objective at each
innerloop iteration. The convergence of the algorithm to a first-order optimal
point is guaranteed by the convergence property of the projected gradient with
Armijo steps [41]. The solution update is computed in line 16, retracted back
onto the constraint manifold (line 17), and the gradient ascent continues
until convergence.

Computational Complexity: The worst-case complexity of Algorithm 2 is
$\mathrm{O}(n^{4})$. This is because the gradient computation (line 19)
involves matrix-vector multiplications, which have $\mathrm{O}(n^{2})$
complexity, and profiling shows this is where most of the time is spent. The
number of backtracking iterations (line 15) and gradient ascent iterations
(line 13) can vary depending on the parameters and the data matrix, but it is
linear in problem size ($\mathrm{O}(n)$) for quadratic objective
$F(u)=u^{\top}M_{d}u$. Lastly, the number of outerloop iterations (line 7)
depends on $\Delta d$ increments (line 26) required to reach $d\geq n$ (at
which point the solution is guaranteed to converge to the binary state). This
is also linear in problem size.

###  IV-C CLIPPER+ Maximal Clique Algorithm

Both Algorithms 1 and 2 are algorithms for finding maximal cliques. The greedy
approach of Algorithm 1 runs fast but gives relatively less accurate estimates
of the maximum clique size when the graph is not sparse. In contrast, the
optimization approach of Algorithm 2 is relatively slower but more accurate.
The main motivation of the proposed CLIPPER+ algorithm is to combine these two
algorithms and thereby combine their relative benefits.

The enable this combination, recall that the core number of a vertex is the
largest integer $k$ such that the degree of the vertex remains non-zero when
all vertices of degree less than $k$ are removed. If a graph contains a clique
of size $k$, then each vertex in the clique must have a degree of $k-1$ or
larger, and therefore a core number of $k-1$ or larger. For this reason, if we
hope to find a larger clique, say of size $k+1$, then only vertices that have
a core number $k$ or higher can be candidates. Using this observation,
CLIPPER+, detailed in Algorithm 3, first runs the greedy algorithm and obtains
a maximal clique (line 4). Assuming this clique has size $k$, if the graph
contains a larger clique, then vertices must have core numbers of $k$ or
larger. Therefore, the algorithm prunes the graph by removing vertices with
core numbers strictly less than $k$ (line 6). The pruning effectively limits
the search space for the optimization (line 10) by reducing the number of
vertices, which improves the runtime.

If the clique recovered by the greedy algorithm is the maximum clique, then
pruning the graph removes all the vertices and therefore we can terminate
early (line 7). This early termination particularly occurs when the graph is
sparse, which leads to a significant speed-up over running the optimization-
based algorithm on the original graph.

When the optimization-based algorithm is required, the initial guess used for
the optimization can be chosen strategically to improve the chance of finding
the maximum clique. This can be done by choosing an initial guess vector that
is a binary complement of the clique found from the greedy approach (line 9),
and thus help the algorithm to converge to a different clique solution.
Lastly, the best solution is selected as the largest clique (line 11).

Computational Complexity: The worst-cast complexity of Algorithm 3 is
$O(n^{4})$. This is because it sequentially combines Algorithms 1 and 2, and
thus inherits the highest worst-case complexity of its components (all other
steps in Algorithm 3 have lower complexity). This basic worst-case complexity
analysis is independent of any input graph structure. The numerical results in
the experimental section provide a good indication of how the complexity
translates to runtimes for various graph sizes in the practical applications.

##  V Experimental Evaluations

Figure 2: The effect of outliers on the consistency graph. (Left) Putative
associations (red: outliers; green: inliers). (Right) Resulting consistency
graph. The higher the outlier percentage, the sparser the graph.

We evaluate the performance of CLIPPER+ (Algorithm 3) in terms of the maximum
clique estimation accuracy and runtime. In addition, we provide ablation
studies of CLIPPER+ by reporting the results of its standalone greedy
(Algorithm 1) and optimization (Algorithm 2) components. We show that CLIPPER+
achieves a performance superior to these standalone components through
combining them.

Algorithms: We test both classical and state-of-the-art maximum clique
estimation algorithms. Classical works include Pelillo [35] and Ding et al.
[36] based on relaxations of the Motzkin-Straus formulation [37]. State-of-
the-art include our implementation of the greedy parallel maximum clique (PMC)
algorithm [2] in Algorithm 1 (which is theoretically equivalent to ROBIN
[19]), the algorithm by Belachew and Gillis [38], and our prior CLIPPER
algorithm [1].

Benchmarks: Our evaluations examine finding the maximum clique on general
graphs, as well as graphs that result from the graph formulation of
synthetic/real-world point cloud registration problems. While general graphs
can have any structure, the registration graphs have certain patterns (e.g.,
sparsity) that significantly affect the results. Fortunately, as we will see,
finding the maximum clique is empirically easier on graphs that result from
registration.

Platform and implementations: All benchmarks are run on a machine with an
Intel Core i9 Processor and 32 GB RAM. The algorithms by Pelillo, Ding, and
Belachew are implemented in Matlab by their original authors. Our prior work
CLIPPER is also implemented in Matlab. The “Greedy” (Algorithm 1), “Optim”
(Algorithm 2), and CLIPPER+ algorithms are implemented in C++, which interface
to Matlab via binary MEX binders for benchmarking.

###  V-A Maximum Clique Benchmark—DIMACS

Figure 3: The maximum clique estimation accuracy ratio on the Stanford Bunny
benchmark (the closer to $1$, the better). Each point corresponds to the
accuracy ratio from a single trial. CLIPPER+ outperforms all algorithms.

Dataset: To evaluate the algorithms on general graphs, we use the DIMACS
benchmark [40], which was introduced in 1996 and has been used widely since
then to benchmark the maximum clique algorithms. The DIMACS dataset consists
of graphs for which finding the maximum clique is challenging. In the interest
of space, we present the result on a subset of smaller graphs in this dataset,
shown in Table I. While DIMACS graphs are relatively small (100-4000
vertices), they are dense and contain more edges compared to graphs resulting
from registration. Despite the dataset being around for more than 25 years,
the maximum clique of some graphs is still unknown, demonstrating the
difficulty of the problem. For instance, on the `C250.9` graph, the (multi-
threaded) PMC exact algorithm [2], used in this work for ground truth
generation and in TEASER [7] for certifiable registration, took $28$ minutes
on our machine to find the maximum clique.

Evaluation: Table I compares the accuracy and runtime of algorithms. The graph
sparsity is defined as $s\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}1-\frac{|E|}{|E_{\max}|}\in[0,1]$, where $|E|$ is the number of
graph edges, and $|E_{\max}|$ is the maximum possible number of edges. The
small values of $s$ show that DIMACS graphs are generally dense. The accuracy
ratio is measured by $r\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny
def}}}{=}}\hat{\omega}/\omega_{\mathrm{gt}}$, where $\hat{\omega}$ is the
clique size found by the algorithm, and $\omega_{\mathrm{gt}}$ is the ground
truth maximum clique size. The ratio of $1$ indicates that the maximum clique
is found, hence, the closer $r$ is to $1$, the more accurate an algorithm is.
CLIPPER+ outperforms all algorithms in the overall accuracy. Unsurprisingly,
the greedy algorithm has the best overall runtime due to its low computational
complexity; however, it has a lower accuracy. By combining the greedy and
optimization algorithms, CLIPPER+ obtains high accuracy and better overall
runtime compared to the optimization-only approach (e.g., a $\sim$10x
reduction of runtime from $29.4$ to $3.6$ milliseconds on `brock200_2`).

###  V-B Registration Benchmark—Stanford Bunny

Figure 4: Mean maximum clique estimation accuracy across different outlier
percentages on the Stanford Bunny benchmark (the closer to $1$, the better).
Figure 5: Mean runtime across different outlier percentages on the Stanford
Bunny benchmark (the lower, the better). Figure 6: A registration instance
between two scans of the 7-Scenes dataset [42]. From left to right: $688$
putative FPFH associations having $93.6\%$ outliers (red lines); corresponding
consistency graph used as input to CLIPPER+ (the maximum clique is
highlighted); inlier associations corresponding to the maximum clique solution
found by CLIPPER+; aligned point clouds using the rotation/translation
computed from the Arun’s method [43] using inliers. TABLE II:  Comparisons of
maximum clique estimation accuracy and runtime of CLIPPER+ on real-world point
cloud registration datasets (bold is best).

Benchmark |  |  |  |  |  |  | Pelillo |  | Ding |  | Belachew |  | CLIPPER |  | Greedy |  | Optim |  | CLIPPER+
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Name | Sequence | $N$ | $\bar{op}$ | $\bar{n}$ | $\bar{s}$ |  | $\bar{t}\downarrow$ | $\bar{r}\uparrow$ |  | $\bar{t}\downarrow$ | $\bar{r}\uparrow$ |  | $\bar{t}\downarrow$ | $\bar{r}\uparrow$ |  | $\bar{t}\downarrow$ | $\bar{r}\uparrow$ |  | $\bar{t}\downarrow$ | $\bar{r}\uparrow$ |  | $\bar{t}\downarrow$ | $\bar{r}\uparrow$ |  | $\bar{t}\downarrow$ | $\bar{r}\uparrow$
7-Scenes | `redkitchen` | 244 | 92.8% | 680 | 0.88 |  | 118.5 | 0.96 |  | 28.3 | 0.9 |  | 73.7 | 0.99 |  | 62.7 | 0.97 |  | 10.1 | 0.89 |  | 119.9 | 0.99 |  | 53.1 | 0.99
Sun3D | `home_at_scan1` | 81 | 88% | 670 | 0.83 |  | 133.1 | 0.97 |  | 17.5 | 0.87 |  | 64.0 | 0.98 |  | 58.2 | 0.98 |  | 10.3 | 0.89 |  | 76.0 | 0.99 |  | 46.4 | 0.99
`hotel_uc_scan3` | 117 | 93.5% | 843 | 0.89 |  | 215.3 | 0.96 |  | 28.9 | 0.89 |  | 116.5 | 0.99 |  | 109.3 | 0.98 |  | 12.3 | 0.89 |  | 159.4 | 0.99 |  | 56.1 | 0.99
`mit_76_1studyroom2` | 96 | 94.2% | 936 | 0.9 |  | 192.9 | 0.97 |  | 32.2 | 0.87 |  | 120.6 | 0.99 |  | 102.8 | 0.97 |  | 14 | 0.88 |  | 180.8 | 0.99 |  | 83.9 | 0.99
ETH | `gazebo_summer` | 116 | 98.3% | 4320 | 0.97 |  | 4255.2 | 0.95 |  | 1275.6 | 0.8 |  | 8467.8 | 0.98 |  | 7978.5 | 0.97 |  | 351.5 | 0.81 |  | 7595.7 | 0.99 |  | 3769 | 0.99
`wood_autumn` | 87 | 99.4% | 5703 | 0.98 |  | 4866.2 | 0.94 |  | 2990.9 | 0.89 |  | 14464 | 0.98 |  | 12774 | 0.97 |  | 564.5 | 0.75 |  | 14758 | 0.99 |  | 14594 | 0.99

$N$: Total number of overlapping point cloud scans on which registration is
performed. $\bar{op}$: Mean of outlier percentages in putative associations.
$\bar{n}$: Mean graph size. $\bar{s}$: Mean graph sparsity. $\bar{t}$: Mean
runtime (milliseconds); the lower, the better. $\bar{r}$: Mean maximum-clique
accuracy ratio ($\hat{\omega}/\omega_{\mathrm{gt}}$); the closer to $1$, the
better.

Using the Stanford Bunny point cloud [44], shown in Fig. 2, we evaluate
CLIPPER+ as an algorithm for global point cloud registration in various
outlier regimes.

Dataset: The (downsampled) Bunny point cloud consists of $1000$ points that
fit in a cube of size $0.2\,\mathrm{m}$. To generate a second point cloud, we
add uniform noise in the range $[-\epsilon/2,\epsilon/2]$ to all points, where
$\epsilon$ is set as the mean distance of all points to their nearest
neighbors in the Bunny point cloud. Additionally, $1000$ outlier points
randomly drawn from a sphere of radius $1\,\mathrm{m}$, centered at the bunny
point cloud, are added to simulate clutter. From the set of all possible
associations between the points in the two point clouds, $200$ associations
are randomly selected from the set of inlier and outlier associations (we keep
this number small to be able to find the ground truth maximum clique). We
consider different outlier ratios (ranging from $0\%$ to $98\%$ in $2\%$
increments), to test scenarios with various data association accuracy, as
shown in Fig. 2. We use the noise $\epsilon$ as the threshold to generate the
consistency graph (according to Section III). The ground truth maximum clique
is found by running the exact PMC algorithm in [2]. Evaluations of maximum
clique estimation accuracy and runtime are performed across $50$ Monte Carlo
runs/graphs for each increment of the outlier percentage.

Evaluation: Fig. 3 compares the accuracy ratio $r$ of algorithms across all
runs and all outlier ratios. CLIPPER+ clearly demonstrates the highest
accuracy because the distribution of $r$ is closest to $1$. Interestingly,
while the greedy and optimization algorithms returned low-accuracy solutions
in some instances, by combining these algorithms CLIPPER+ obtains an accuracy
beyond these standalone components.

Fig. 4 shows the mean of the accuracy ratio $r$ (averaged across $50$ Monte
Carlo runs at each outlier percentage increment) versus the outlier
percentage. The accuracy of the greedy approach is low in low-outlier regimes
(which have dense consistency graphs, see Fig. 2). As the outlier percentage
increases and the graph becomes sparser, the accuracy of the greedy method
improves. Both CLIPPER+ and optimization algorithms retain a high accuracy
ratio across all outlier percentages.

The mean runtime of the algorithms is shown in Fig. 5. The runtime of the
optimization method remains roughly the same, while the runtime of the greedy
method improves as the outlier ratio grows and the graph becomes sparser. The
runtime of CLIPPER+ in low outlier regimes is roughly equal to the compound
runtimes of its greedy and optimization components because in such regimes
pruning the graph based on core numbers does not remove a significant number
of vertices (if any). However, as the outlier ratio increases and the graph
becomes sparser, pruning removes more vertices and its speed-up effect becomes
more apparent. This can be seen around the $20\%$ outlier ratio, where
CLIPPER+ becomes faster than the optimization method, and its speed improves
further as the outlier percentage increases.

###  V-C Registration Benchmark—Real-World Point Clouds

Datasets: We use sequences in the real-world 7-Scenes [42], Sun3D [45], and
ETH [46] datasets (similar to 3DMatch [15] and 3DSmoothNet [47] evaluations).
Sun3D and 7-Scenes are dense indoor RGB-D point clouds, while ETH is an
outdoor LiDAR point cloud. In each sequence, we consider pairs of point clouds
(or scans) that have an overlap. To increase registration speed, we downsample
the point clouds by discretizing the 3D space into cubes of size
$\epsilon=0.05\,\mathrm{m}$ for 7-Scenes and Sun3D datasets, and
$\epsilon=0.1\,\mathrm{m}$ for the ETH dataset, and using the mean of the
point coordinates in each cube as a single-point representative. For each
downsampled point, FPFH descriptor vectors [14] are computed and associated
bilaterally based on their $l_{2}$ norm distance using the k-nearest neighbors
algorithm. To generate the ground truth for our evaluations, we use the exact
maximum clique algorithm of [2] to find the largest set of geometrically
consistent associations in these putative FPFH associations. We store results
if this maximum clique solution correctly registers the point clouds according
to the ground truth provided by the datasets—the maximum clique/likelihood
solution may register points wrongly due to practical limitations such as
repetitive patterns (perceptual aliasing), insignificant overlap between the
point clouds, and lack of any inlier associations caused by downsampling and
FPFH inaccuracies. For evaluations, the consistency graph is generated
according to Section III from putative FPFH associations, using the
downsampling $\epsilon$ as the consistency threshold. An evaluation instance
is shown in Fig. 6.

Evaluation: Table II presents the evaluation results for all datasets and
algorithms. CLIPPER+ outperforms all algorithms in accuracy on all
datasets/sequences. The greedy algorithm has the smallest runtime at the
expense of the lowest overall accuracy. The standalone optimization algorithm
has similar accuracy to CLIPPER+. However, except on the last sequence, it is
around $2$x slower. This demonstrates the advantage of CLIPPER+ over its
standalone greedy and optimization components.

Lastly, we point out the high outlier ratios of FPFH associations (e.g., on
average, $99.4\%$ of associations in the `wood_autumn` sequence are outliers).
Due to the maximum clique solution correctly registering the point clouds in
our benchmark, the accuracy ratio $r$ shows the point cloud registration
success rate. Thus, CLIPPER+ correctly registers the point clouds in $99\%$ of
the trials despite extreme outlier percentages. This is what distinguishes
CLIPPER+ from existing robust registration frameworks (such as RANSAC [16]),
that can fail in these high-outlier regimes.

##  VI Conclusion and Future Directions

We presented CLIPPER+, a maximal-clique-finding algorithm for unweighted
graphs that enables robust global registration in robotics and computer vision
applications. Future work includes investigating alternative optimization
methods such as second-order or quasi-Newton methods, and an extension to the
weighted graphs based on our previous work [1] (designed for weighted graphs)
and an extension of core numbers to weighted settings (studied in [48]). We
also plan to integrate the algorithm in point cloud registration pipelines [7]
as the outlier rejection module to improve the runtime and outlier rejection
capacity.

## References

  * [1] P. C. Lusk, K. Fathian, and J. P. How, “CLIPPER: A graph-theoretic framework for robust data association,” in _IEEE Int. Conf. on Robotics and Automation_ , 2021, pp. 13 828–13 834.
  * [2] R. A. Rossi, D. F. Gleich, and A. H. Gebremedhin, “Parallel maximum clique algorithms with applications to network analysis,” _SIAM Journal on Scientific Computing_ , vol. 37, no. 5, pp. C589–C616, 2015.
  * [3] P. Y. Lajoie, B. Ramtoula, Y. Chang, L. Carlone, and G. Beltrame, “DOOR-SLAM: Distributed, online, and outlier resilient slam for robotic teams,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 1656–1663, 2020.
  * [4] J. G. Mangelson, D. Dominic, R. M. Eustice, and R. Vasudevan, “Pairwise consistent measurement set maximization for robust multi-robot map merging,” in _IEEE Int. Conf. on Robotics and Automation_ , 2018, pp. 2916–2923.
  * [5] J. Ankenbauer, K. Fathian, and J. P. How, “View-invariant localization using semantic objects in changing environments,” _arXiv preprint arXiv:2209.14426_ , 2022.
  * [6] A. P. Bustos and T.-J. Chin, “Guaranteed outlier removal for point cloud registration with correspondences,” _IEEE trans. on pattern analysis & machine intelligence_, vol. 40, no. 12, pp. 2868–2882, 2017.
  * [7] H. Yang, J. Shi, and L. Carlone, “TEASER: Fast and certifiable point cloud registration,” _IEEE Transactions on Robotics_ , vol. 37, no. 2, pp. 314–333, 2020.
  * [8] Y. Li, L. Gu, and T. Kanade, “Robustly aligning a shape model and its application to car alignment of unknown pose,” _IEEE trans. on pattern analysis & machine intelligence_, vol. 33, no. 9, pp. 1860–1876, 2011.
  * [9] C. R. Qi, H. Su, K. Mo, and L. J. Guibas, “PointNet: Deep learning on point sets for 3d classification and segmentation,” in _IEEE conf. on comp. vision & pattern recognition_, 2017, pp. 652–660.
  * [10] K. Fathian, K. Khosoussi, Y. Tian, P. Lusk, and J. P. How, “CLEAR: A consistent lifting, embedding, and alignment rectification algorithm for multiview data association,” _IEEE Transactions on Robotics_ , vol. 36, no. 6, pp. 1686–1703, 2020.
  * [11] P. C. Lusk, K. Fathian, and J. P. How, “MIXER: Multiattribute, multiway fusion of uncertain pairwise affinities,” _IEEE Robotics and Automation Letters_ , vol. 8, no. 5, pp. 2462–2469, 2023.
  * [12] L. Rakai, H. Song, S. Sun, W. Zhang, and Y. Yang, “Data association in multiple object tracking: A survey of recent techniques,” _Expert Systems with Applications_ , vol. 192, p. 116300, 2022.
  * [13] P. J. Besl and N. D. McKay, “Method for registration of 3-D shapes,” in _Sensor fusion IV: control paradigms and data structures_ , vol. 1611\. Int. Society for Optics and Photonics, 1992, pp. 586–606.
  * [14] R. B. Rusu, N. Blodow, and M. Beetz, “Fast point feature histograms (FPFH) for 3D registration,” in _IEEE ICRA_ , 2009, pp. 3212–3217.
  * [15] A. Zeng, S. Song, M. Nießner, M. Fisher, J. Xiao, and T. Funkhouser, “3DMatch: Learning Local Geometric Descriptors from RGB-D Reconstructions,” in _IEEE conf. on comp. vision and pattern recognition_ , 2017, pp. 1802–1811.
  * [16] M. A. Fischler and R. C. Bolles, “Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography,” _Communications of the ACM_ , vol. 24, no. 6, pp. 381–395, 1981.
  * [17] H. Yang, P. Antonante, V. Tzoumas, and L. Carlone, “Graduated non-convexity for robust spatial perception: From non-minimal solvers to global outlier rejection,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 1127–1134, 2020.
  * [18] P. C. Lusk, D. Parikh, and J. P. How, “GraffMatch: Global Matching of 3D Lines and Planes for Wide Baseline LiDAR Registration,” _IEEE Robotics and Automation Letters_ , vol. 8, no. 2, pp. 632–639, 2022.
  * [19] J. Shi, H. Yang, and L. Carlone, “ROBIN: a graph-theoretic approach to reject outliers in robust estimation using invariants,” in _IEEE Int. Conf. on Robotics and Automation_ , 2021, pp. 13 820–13 827.
  * [20] G. Izatt, H. Dai, and R. Tedrake, “Globally optimal object pose estimation in point clouds with mixed-integer programming,” in _Robotics Research: The 18th Int. Symposium_. Springer, 2020, pp. 695–710.
  * [21] J. Yang, H. Li, D. Campbell, and Y. Jia, “Go-icp: A globally optimal solution to 3d icp point-set registration,” _IEEE Trans. on pattern analysis & machine intelligence_, vol. 38, no. 11, pp. 2241–2254, 2015.
  * [22] T. J. Chin and D. Suter, “The maximum consensus problem: recent algorithmic advances,” _Synthesis Lectures on Comp. Vision_ , vol. 7, no. 2, pp. 1–194, 2017.
  * [23] D. Barath and J. Matas, “Graph-cut ransac,” in _& IEEE conf. on comp. vision and pattern recognition_, 2018, pp. 6733–6741.
  * [24] M. J. Black and A. Rangarajan, “On the unification of line processes, outlier rejection, and robust statistics with applications in early vision,” _Int. journal of comp. vision_ , vol. 19, no. 1, pp. 57–91, 1996.
  * [25] A. Blake and A. Zisserman, _Visual reconstruction_. MIT press, 1987.
  * [26] P. Agarwal, G. D. Tipaldi, L. Spinello, C. Stachniss, and W. Burgard, “Robust map optimization using dynamic covariance scaling,” in _IEEE Int. Conf. on Robotics and Automation_ , 2013, pp. 62–69.
  * [27] J. L. Schonberger and J.-M. Frahm, “Structure-from-motion revisited,” in _IEEE conf. on comp. vision & pattern recognition_, 2016, pp. 4104–4113.
  * [28] R. Raguram, J.-M. Frahm, and M. Pollefeys, “A comparative analysis of ransac techniques leading to adaptive real-time random sample consensus,” in _European Conf. on Comp. Vision_. Springer, 2008, pp. 500–513.
  * [29] T. Bailey, E. M. Nebot, J. Rosenblatt, and H. F. Durrant-Whyte, “Data association for mobile robot navigation: A graph theoretic approach,” in _IEEE Int. Conf. on Robotics and Automation_ , vol. 3, 2000, pp. 2512–2517.
  * [30] M. Leordeanu and M. Hebert, “A spectral technique for correspondence problems using pairwise constraints,” in _IEEE Int. Conf. on Comp. Vision_ , vol. 2, 2005, pp. 1482–1489.
  * [31] O. Enqvist, K. Josephson, and F. Kahl, “Optimal correspondences from pairwise constraints,” in _IEEE Int. Conf. on Comp. Vision_ , 2009, pp. 1295–1302.
  * [32] A. P. Bustos, T.-J. Chin, F. Neumann, T. Friedrich, and M. Katzmann, “A practical maximum clique algorithm for matching with pairwise constraints,” _arXiv preprint arXiv:1902.01534_ , 2019.
  * [33] X. Zhang, J. Yang, S. Zhang, and Y. Zhang, “3D Registration with Maximal Cliques,” in _IEEE/CVF Conf. on Comp. Vision and Pattern Recognition_ , 2023, pp. 17 745–17 754.
  * [34] Q. Wu and J.-K. Hao, “A review on algorithms for maximum clique problems,” _European Journal of Operational Research_ , vol. 242, no. 3, pp. 693–709, 2015.
  * [35] M. Pelillo, “Relaxation labeling networks for the maximum clique problem,” _Journal of Artificial Neural Networks_ , vol. 2, no. 4, pp. 313–328, 1995\.
  * [36] C. Ding, T. Li, and M. I. Jordan, “Nonnegative matrix factorization for combinatorial optimization: Spectral clustering, graph matching, and clique finding,” in _IEEE Int. Conf. on Data Mining_ , 2008, pp. 183–192.
  * [37] T. S. Motzkin and E. G. Straus, “Maxima for graphs and a new proof of a theorem of Turan,” _Canadian Journal of Mathematics_ , vol. 17, pp. 533–540, 1965.
  * [38] M. T. Belachew and N. Gillis, “Solving the maximum clique problem with symmetric rank-one non-negative matrix approximation,” _Journal of Optimization Theory & Applications_, vol. 173, no. 1, pp. 279–296, 2017.
  * [39] V. Batagelj and M. Zaversnik, “An O(m) algorithm for cores decomposition of networks,” _arXiv preprint cs/0310049_ , 2003.
  * [40] D. S. Johnson and M. A. Trick, _Cliques, coloring, and satisfiability: second DIMACS implementation challenge, October 11-13, 1993_. American Mathematical Soc., 1996, vol. 26.
  * [41] D. P. Bertsekas, “Nonlinear programming,” _Journal of the Operational Research Society_ , vol. 48, no. 3, pp. 334–334, 1997.
  * [42] J. Shotton, B. Glocker, C. Zach, S. Izadi, A. Criminisi, and A. Fitzgibbon, “Scene coordinate regression forests for camera relocalization in rgb-d images,” in _IEEE conf. on comp. vision and pattern recognition_ , 2013, pp. 2930–2937.
  * [43] K. S. Arun, T. S. Huang, and S. D. Blostein, “Least-squares fitting of two 3-d point sets,” _IEEE Transactions on pattern analysis and machine intelligence_ , no. 5, pp. 698–700, 1987.
  * [44] B. Curless and M. Levoy, “A volumetric method for building complex models from range images,” in _23rd annual conf. on comp. graphics and interactive techniques_ , 1996, pp. 303–312.
  * [45] J. Xiao, A. Owens, and A. Torralba, “SUN3D: A database of big spaces reconstructed using SfM and object labels,” in _IEEE Int. conf. on comp. vision_ , 2013, pp. 1625–1632.
  * [46] F. Pomerleau, M. Liu, F. Colas, and R. Siegwart, “Challenging data sets for point cloud registration algorithms,” _The Int. Journal of Robotics Research_ , vol. 31, no. 14, pp. 1705–1711, 2012.
  * [47] Z. Gojcic, C. Zhou, J. D. Wegner, and W. Andreas, “The perfect match: 3d point cloud matching with smoothed densities,” in _Int. conf. on comp. vision and pattern recognition_ , 2019.
  * [48] B. Liu and F. Zhang, “Incremental algorithms of the core maintenance problem on edge-weighted graphs,” _IEEE Access_ , vol. 8, pp. 63 872–63 884, 2020\.