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

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

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

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

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

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

\textbf{Contributions:} We present the following contributions.
\begin{enumerate}
	\item To the best of our knowledge, our proposed approach is a first attempt to formulate
	\emph{higher-order matching} between two sets of points, given partial or incomplete information, as a matching between two \emph{random clique complexes}. We also propose an efficient matching algorithm and study both its time and storage complexity.

	\item (i) We provide new bounds on the concentration inequality of eigenvalues of the QAP trace formulation for random symmetric matrices, (ii) we give tighter concentration inequality bounds on the largest eigenvalue for the Lawler QAP formulation on random matrices, in the context of \emph{affinity matrices} that are used by some earlier works.  Furthermore, we theoretically analyze and discuss the robustness of affinity-matrix based schemes to missing points, and
	(iii) we perform asymptotic analysis on the worst to best case ratio of a QAP solution for our higher-dimensional clique adjacency matrices in the \emph{clique percolation regime}~\cite{Bollobas2009},
	where the entries follow a \emph{Poisson} distribution.

	\item Finally, we present a comprehensive empirical study that compares our method's matching accuracy to that of a diverse set of matching approaches~\cite{ZhouD16,zhou2013deformable,cho2010reweighted,feizi2016spectral,Leordeanu2005,cour2007balanced,pachauri2013solving,gold1996graduated,Kuhn1955,leordeanu2009integer,zass2008probabilistic,li2013object,duchenne2011tensor}.
	We conducted our experiments on both synthetic and well-known hard real-world datasets that span across affine/non-affine transformations, severe occlusions, and clutter. Our study reveals much better accuracy for the popular datasets against several of the state-of-the-art matching methods.
\end{enumerate}

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

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

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

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

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

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

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

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

\begin{algorithm}[tbp]
	\caption{Matching Random Clique Complexes}
	\label{alg:matching}
	\begin{flushleft}
		\textbf{Input:} $\mathcal{X}(G) = \{  G^{(k,l)}   \}_{k=0}^h$ and
		$\mathcal{X}(G') = \{ G'^{(k,l)}   \}_{k=0}^h$
	\end{flushleft}
	\begin{algorithmic}[1]

		\FOR{ $k = h \dots 0$}
		\STATE Let $M, M'$ be the total number of $(k+1)$-cliques in $G^{(k,l)}$
		and $G'^{(k,l)}$, respectively

		\STATE $\mathcal{L} := \{   c_i^{(k)}  \}_{i=0}^{M-1} $  \COMMENT {list of barycenters}

		\FOR{ $i = 0 \dots M-1$}
		\STATE $\mathcal{N}_i := \mathcal{N}_i \cup \left\{  g^{(k,l)}_{(x,:)} \mid x=i,   g^{(k,l)}_{(x,y)} \neq 0   \right\}  $
		\STATE $\mathcal{N} := \mathcal{N} \cup \{  \mathcal{N}_i \} $ \COMMENT {clique neighborhoods}
		\ENDFOR

		\FOR{ $i = 0 \dots M-1$}
		\STATE $\alpha_i := [ \alpha_1, \dots, \alpha_{M-1}]^T$
		\STATE $\alpha := \alpha \cup \{ \alpha_i \}$  \COMMENT{affine weight vectors}
		\ENDFOR

		\STATE Repeat steps $3$--$11$ on $G'^{(k,l)}$ for $\mathcal{L}',\mathcal{N}'$ and $\alpha'$.

		\STATE Build cost matrix $C^{(k)}$ from weights vectors $\alpha, \alpha'$
		\STATE $X^{*}_k := $ \textbf{Kuhn-Munkres }($G^{(k,l)},G'^{(k,l)},C^{(k)} $)
		\ENDFOR
	\end{algorithmic}
	\begin{flushleft}
		\textbf{Return:}  $ \{ X^{*}_0, \dots, X^{*}_h \}$ \text{ \# set of permutation matrices}
	\end{flushleft}
\end{algorithm}

\subsection{Problem Setup}
\label{subsec:setup}
The problem of matching random clique complexes each of dimension $h$, is the estimation of a set of
optimal bijective maps of the form $\mathcal{M}_i: \mathcal{X}^{(i)}(G) \rightarrow \mathcal{X}^{(i)}(G')$, for all $i \leq h$, subject to assignment constraints.
This can be formulated as a \emph{constrained quadratic assignment problem}, which can later be relaxed to a linear programming optimization problem.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

	\end{enumerate}
\end{theorem}

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

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

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

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

	As $\vertiii{X}$ is non-negative, we get $|\lambda_i| \leq \vertiii{A}$.
	Thus, every eigenvalue of $A$ is upper bounded by the matrix norm $\vertiii{A}$.
	Applying the triangle inequality, we get that $\delta_i(A) = |\lambda_{i+1}(A) - \lambda_i(A)| \leq 2 \vertiii{A}$, which completes the proof.
\end{proof}

\begin{proof}
	For ease of understanding, we drop the $(A)$ as it is obvious from context. Let $\lambda_i$ be the $i$-th eigenvalue of $A$, and let $x_i \neq 0$ be its corresponding eigenvector. From $A x_i = \lambda_i x_i$, we have
	\begin{align*}
	A X_i = \lambda_i X_i, \text{\hspace{3em}where }
	X_i:= \begin{bmatrix}[c|c|c]
	x_i & \dots & x_i
	\end{bmatrix}
	\in \mathcal{M}_n \setminus \{0\}
	\end{align*}
	It follows,
	\begin{align*}
	|\lambda_i|  \vertiii{X_i} \vertiii{X_i}  = \vertiii{\lambda_i X_i}
	=  \vertiii{AX_i} \leq \vertiii{A} \vertiii{X_i} .
	\end{align*}

	As $\vertiii{X}$ is non-negative, we get $|\lambda_i| \leq \vertiii{A}$.
	Thus, every eigenvalue of $A$ is upper bounded by the matrix norm $\vertiii{A}$.
	Applying the triangle inequality, we get that $\delta_i(A) = |\lambda_{i+1}(A) - \lambda_i(A)| \leq 2 \vertiii{A}$, which completes the proof.
\end{proof}
\fi

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

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

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

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

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

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

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

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

\begin{figure*}[t!]
	\centering
	\subfigure[]{%
		\label{fig:zero11}%
		\includegraphics[width=27mm]{Experiments/Transformation/house0.png}}%
	\qquad
	\subfigure[]{%
		\label{fig:first11}%
		\includegraphics[width=27mm]{Experiments/Transformation/HouseRotate180.png}}%
	\qquad
	\subfigure[]{%
		\label{fig:second11}%
		\includegraphics[width=27mm]{Experiments/Transformation/HouseMirrorx_marked.png}}%
	\qquad
	\subfigure[]{%
		\label{fig:third11}%
		\includegraphics[width=27mm,height=27mm]{Experiments/Transformation/HouseScale05_15.png}}%
	\qquad
	\subfigure[]{%
		\label{fig:fourth11}%
		\includegraphics[width=27mm,height=27mm]{Experiments/Transformation/HouseShear05_marked.png}}%
	\vspace{-5mm}
	\caption{(a) Original House Frame, (b)--(e) Four transformations on house frame: Rotation, Reflection, Scaling and Shear (green markers show true matching case with the original frame (a)).}
	\label{fig:HouseTransform}
\end{figure*}
\begin{table*}[h]
	{\tiny
		\hfill{}
		\begin{tabular}{l||l|l|l|l|l|l|l|l|l|l}
			\hline
			\multicolumn{1}{l||}{\textbf{Algorithms}} & \multicolumn{2}{|l|}{\textbf{$20^{\circ}$} Rotation} & \multicolumn{2}{|l|}{\textbf{$60^{\circ}$} Rotation}&  \multicolumn{2}{|l|}{\textbf{Reflection}} & \multicolumn{2}{|l|}{\textbf{Scaling}} & \multicolumn{2}{|l}{\textbf{Shear}}\\

			\cline{1-11}
			& 20\% & 40\% & 20\% & 40\% & 20\% & 40\% & 20\% & 40\% & 20\% & 40\% \\
			\cline{2-11}
			\hline
			\hline
			OurMethod & \textbf{0.01 $\pm$ 0.0} & \textbf{0.01 $\pm$ 0.0} & \textbf{0.05 $\pm$ 0.0} & \textbf{0.03 $\pm$ 0.0} & \textbf{0.0 $\pm$ 0.0} & \textbf{0.0 $\pm$ 0.0} & \textbf{0.1 $\pm$ 0.0} & \textbf{0.2 $\pm$ 0.1} & \textbf{0.0 $\pm$ 0.0} & \textbf{0.0 $\pm$ 0.0}\\
			EigenAlign & 63.97 $\pm$ 1.0 & 70.19 $\pm$ 1.0 & 65.39 $\pm$ 1.5 & 70.88 $\pm$ 1.9 & 62.5 $\pm$ 0.4 & 64.97 $\pm$ 0.2 & 62.37 $\pm$ 0.6 & 65.42 $\pm$ 1.1 & 61.04 $\pm$ 0.2 & 62.36 $\pm$ 0.7\\

			FGM & 3.6 $\pm$ 0.5 & 7.4 $\pm$ 0.5 & 18.0 $\pm$ 0.0 & 36.4 $\pm$ 0.5 & \textbf{0.0 $\pm$ 0.0} & \textbf{0.0 $\pm$ 0.0} & 2.20 $\pm$ 1.0 & 3.4 $\pm$ 0.5 & \textbf{0.0 $\pm$ 0.0} & \textbf{0.0 $\pm$ 0.0}\\
			LAI-LP & 40.90 $\pm$ 0.9 & 43.07 $\pm$ 1.2 & 49.9 $\pm$ 1.2 & 61.16 $\pm$ 0.7 & 49.07 $\pm$ 0.9 & 59.41 $\pm$ 0.8 & 43.12 $\pm$ 5.5 & 45.91 $\pm$ 3.2 & 39.04 $\pm$ 0.7 & 38.27 $\pm$ 0.6\\
			PermSync & 13.36 $\pm$ 1.1 & 16.48 $\pm$ 0.4 & 25.48 $\pm$ 1.4 & 41.59 $\pm$ 0.7 & 26.24 $\pm$ 3.4 & 42.31 $\pm$ 5.1 & 12.15 $\pm$ 0.7 & 13.7 $\pm$ 0.6 & 10.19 $\pm$ 0.8 & 6.23 $\pm$ 3.2\\

			RRWM & 2.0 $\pm$ 0.0 & 4.0 $\pm$ 0.0 & 14.0 $\pm$ 0.0 & 27.0 $\pm$ 0.0 & \textbf{0.0 $\pm$ 0.0} & \textbf{0.0 $\pm$ 0.0} & 5.0 $\pm$ 2.0 & 10.0 $\pm$ 1.5 & \textbf{0.0 $\pm$ 0.0} & \textbf{0.0 $\pm$ 0.0}\\
			Tensor & 6.88 $\pm$ 0.0 & 14.39 $\pm$ 0.8 & 19.24 $\pm$ 0.6 & 37.95 $\pm$ 0.8 & 17.1 $\pm$ 0.5 & 34.69 $\pm$ 0.4 & 2.91 $\pm$ 0.6 & 6.19 $\pm$ 0.7 & \textbf{0.0 $\pm$ 0.0} & \textbf{0.0 $\pm$ 0.0}\\
			IPFP & 6.8 $\pm$ 1.5 & 11.6 $\pm$ 0.5 & 18.2 $\pm$ 0.5 & 35.0 $\pm$ 0.0 & 1.0 $\pm$ 0.0 & 1.0 $\pm$ 0.0 & 7.4 $\pm$ 2.0 & 15.6 $\pm$ 1.5 & 0.8 $\pm$ 0.5 & 0.6 $\pm$ 0.5\\
			PM & 43.2 $\pm$ 1.0 & 49.4 $\pm$ 1.5 & 46.4 $\pm$ 1.0 & 56.8 $\pm$ 0.5 & 37.0 $\pm$ 0.0 & 37.0 $\pm$ 0.0 & 45.6 $\pm$ 0.5 & 53.4 $\pm$ 0.5 & 37.8 $\pm$ 0.5 & 38.4 $\pm$ 0.5\\

			SMAC & 12.6 $\pm$ 0.5 & 20.0 $\pm$ 0.0 & 19.0 $\pm$ 0.0 & 33.0 $\pm$ 0.0 & 4.0 $\pm$ 0.0 & 4.0 $\pm$ 0.0 & 13.8 $\pm$ 2.5 & 23.4 $\pm$ 0.5 & 4.0 $\pm$ 0.0 & 4.2 $\pm$ 0.5\\
			SM & 32.6 $\pm$ 0.5 & 39.4 $\pm$ 1.5 & 37.8 $\pm$ 1.0 & 49.6 $\pm$ 1.0 & 25.0 $\pm$ 0.0 & 25.0 $\pm$ 0.0 & 33.0 $\pm$ 1.5 & 40.4 $\pm$ 0.5 & 25.2 $\pm$ 0.5 & 24.8 $\pm$ 0.5\\
			GA & 34.0 $\pm$ 1.5 & 37.0 $\pm$ 1.0 & 38.4 $\pm$ 1.0 & 47.0 $\pm$ 0.0 & 30.0 $\pm$ 0.0 & 30.0 $\pm$ 0.0 & 37.6 $\pm$ 1.0 & 45.4 $\pm$ 0.5 & 29.0 $\pm$ 0.0 & 28.4 $\pm$ 0.5\\
			Munkres & 35.25 $\pm$ 1.5 & 36.47 $\pm$ 0.9 & 44.75 $\pm$ 0.8 & 55.13 $\pm$ 1.3 & 46.05 $\pm$ 1.2 & 57.25 $\pm$ 2.2 & 37.32 $\pm$ 1.5 & 37.81 $\pm$ 1.8 & 32.84 $\pm$ 0.3 & 31.27 $\pm$ 0.8\\
			\hline
	\end{tabular}}
	\hfill{}
	\caption{Error (\%) of transformation on CMU House: inserted $20\%$ and $40\%$ impurity in CMU House frame sequence randomly for rotation ($20^{\circ}$, $60^{\circ}$), reflection, scaling and shear. Minimum error (\%) is shown in bold. Matching is computed for $111$ frames from the $1^{st}$ frame to the other $110$ frames. Our method shows best performance among all the methods.}
	\label{tb:Transformation}
\end{table*}

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

\subsection{Effect of Affine Transformations}
\textbf{Simulated Dataset:} We perform affine transformations on \href{http://vasc.ri.cmu.edu/idb/html/motion/house/index.html}{CMU House}, which is a sequence of $N$ frames extracted from a video. More specifically, we uniformly sample frames (at $20\%$ and $40\%$) and perform affine transformations on the selected frames to distort them. Figure~\ref{fig:HouseTransform} shows examples of affine transformations on house frame sequences. Table \ref{tb:Transformation} shows the comparative error in matching for all the algorithms. We now describe each affine transformation as performance metrics in our experiments.
\begin{center}
	\begin{table*}
		{\tiny
			\hfill{}
			\begin{tabular}{l||l|l|l|l|l|l}
				\hline
				\textbf{Methods} & \textbf{Car} & \textbf{Bike} & \textbf{Butterfly} & \textbf{Magazine} & \textbf{Building} & \textbf{Book}\\
				\hline
				\hline
				OurMethod & \textbf{4.14 $\pm$ 2.45}/ 7.13 & \textbf{3.15 $\pm$ 0.32}/ 6.96 & \textbf{3.89 $\pm$ 0.23}/ 14.76 & \textbf{0.48 $\pm$ 0.02}/ 43.99 & \textbf{4.17 $\pm$ 0.32}/ 12.65 & \textbf{22.20 $\pm$ 1.16}/ 14.86\\

				EigenAlign & 60.68 $\pm$ 0.29/ 19.37 & 57.44 $\pm$ 0.37/ 19.32 & 66.57 $\pm$ 0.0/ 26.13 & 43.23 $\pm$ 0.0/ 93.60 & 90.51 $\pm$ 0.0/ 2.64 & 98.41 $\pm$ 0.0/ 8.29\\
				FGM & 55.51 $\pm$ 0.0/ 1793.9 & 48.17 $\pm$ 0.0/ 2013.7 & 16.12 $\pm$ 0.0/ 674.94 & \textbf{0.0 $\pm$ 0.0}/ 777.55 & 74.87 $\pm$ 0.05/ 2530.5 & 97.54 $\pm$ 0.01/ 4293.9\\
				LAI-LP & 73.06 $\pm$ 0.23/ 152.47 & 42.00 $\pm$ 0.24/ 154.15 & 49.54 $\pm$ 0.0/ 161.06 & 88.73 $\pm$ 0.1/ 184.153 & 87.98 $\pm$ 0.0/ 33.06 & 96.38 $\pm$ 0.0/ 14.71\\
				PermSync & 10.63 $\pm$ 0.0/ 0.45 & 8.90 $\pm$ 0.0/ 0.46 & 46.93 $\pm$ 0.0/ 0.43 & 79.88 $\pm$ 0.0/ 1.08 & 64.00 $\pm$ 0.0/ 0.22 & 70.00 $\pm$ 0.0/ 0.48\\

				RRWM & 60.91 $\pm$ 0.0/ 4.96 & 54.53 $\pm$ 0.0/ 4.83 & 30.99 $\pm$ 0.0/ 8.53 & 1.98 $\pm$ 0.0/ 18.09 & 72.87 $\pm$ 0.01/ 7.98 & 87.04 $\pm$ 0.0/ 21.84\\
				Tensor & 24.37 $\pm$ 0.9/ 93.36 & 15.07 $\pm$ 1.0/ 93.97 & \textbf{1.07 $\pm$ 0.17}/ 107.93 & \textbf{0.0 $\pm$ 0.0}/ 182.07 & 43.24 $\pm$ 2.98/ 40.21 & 32.35 $\pm$ 0.15/ 40.41\\
				IPFP & 65.13 $\pm$ 0.0/ 6.35 & 60.81 $\pm$ 0.0/ 6.28 & 40.90 $\pm$ 0.0/ 8.43 & 3.94 $\pm$ 0.0/ 12.31 & 76.19 $\pm$ 0.0/ 4.65 & 87.74 $\pm$ 0.0/ 8.90\\

				PM & 74.63 $\pm$ 0.0/ 7.07 & 71.93 $\pm$ 0.0/ 4.90 & 70.27 $\pm$ 0.0/ 0.94 & 48.82 $\pm$ 0.0/ 1.69 & 83.79 $\pm$ 0.02/ 2.98 & 91.43 $\pm$ 0.24/ 0.44\\
				SMAC & 70.00 $\pm$ 0.0/ 5.75 & 67.36 $\pm$ 0.0/ 5.52 & 50.53 $\pm$ 0.0/ 4.15 & 5.52 $\pm$ 0.0/ 6.90 & 78.56 $\pm$ 0.22/ 1.94 & 87.88 $\pm$ 0.11/ 3.42 \\
				SM & 68.54 $\pm$ 0.0/ 3.34 & 67.18 $\pm$ 0.0/ 3.47 & 65.96 $\pm$ 0.0/ 3.32 & 34.16 $\pm$ 0.0/ 4.93 & 80.27 $\pm$ 0.07/ 1.88 & 88.66 $\pm$ 0.09/ 2.10\\
				GA & 65.06 $\pm$ 0.0/ 4.53 & 64.60 $\pm$ 0.0/ 4.61 & 61.58 $\pm$ 0.0/ 4.08 & 31.62 $\pm$ 0.0/ 5.83 & 77.20 $\pm$ 0.27/ 3.51 & 87.02 $\pm$ 0.16/ 32.28\\
				Munkres & 33.71 $\pm$ 0.0/ 1.52 & 29.99 $\pm$ 0.0/ 1.49 & 51.87 $\pm$ 0.0/ 1.39 & 79.69 $\pm$ 0.0/ 2.45 & 74.00 $\pm$ 0.0/ 0.73 & 92.00 $\pm$ 0.0/ 1.16\\

				\hline
		\end{tabular}}
		\hfill{}
		\caption{Car, Motorbike~\cite{cho2013learning}, Butterfly, Magazine~\cite{jiang2011linear}, Building and Books error (\%) for pairwise matchings. Computation time (in seconds) is mentioned after the "/" in the above Table.}
		\label{tb:n2}
	\end{table*}
\end{center}
\vspace{-0.8cm}
\textbf{Rotation:}
Figure~\ref{fig:first11} shows a $180^{\circ}$ rotated version of the original house frame (Figure~\ref{fig:zero11}). Table \ref{tb:Transformation} shows errors in matching when
$20$\% of the frames are rotated by both $20^{\circ}$ and $60^{\circ}$, respectively
and when the same transformations are applied to $40$\% of the frames. As the percentage of transformed frames with greater degree increases, we note a substantial increase in error for other methods in comparison to our method's error increase.

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

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

\textbf{Shearing:}
We randomly apply shearing on house in one of the directions with shear factor $0.5$ (shown in Figure~\ref{fig:fourth11}) and measured the performance shown in Table \ref{tb:Transformation}. In addition to our method, we find that affinity-based algorithms also produce robust matchings.

\subsection{Effect of Incomplete and Occluded Landmarks}
To understand the effect of occlusions, we took two real-world datasets, i.e., \emph{Books} and \emph{Building}~\cite{pachauri2013solving} with severe occlusions which are scenes of the same 3D object taken from arbitrary camera angles. These datasets have widely been used in \emph{Structure from Motion (SfM)} problems and are known to be difficult for matching.
Focusing our attention to the last two columns of Table \ref{tb:n2}, it is evident that our method gives the best results.
\begin{figure}[t!]
	\centering
	\subfigure{%

		\includegraphics[trim={0 0 0 0.45cm},clip,width=70mm,height=20mm]{Experiments/Books/BookMatch2_17}}%
	\qquad

	\subfigure{%

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

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

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

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

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

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


\section{Introduction}

Dark matter is one of the persisting problems in modern
theoretical particle physics and cosmology, and up to date, no
experimental verification occurred. Observational data coming from
galactic scale structures, like the bullet cluster, or from
galactic rotation curves, strongly indicate that dark matter seems
to control the dynamics of galaxies rotation and also seems to
control the collision of galaxies. In addition, the observation of
the bullet cluster, strongly indicates that dark matter is a
particle, and halos of dark matter particles accompany the visible
structure of a galaxy. Although, the dark matter effects can be
mimicked by modified gravity itself
\cite{reviews4,reviews5,reviews6,reviews1,reviews2,reviews3}, see
also  \cite{Capozziello:2006ph,Nojiri:2016vhu}, observations seem
to favor the particle nature of dark matter for the moment at
least. Up to present date, the dark matter searches utterly failed
to identify any dark matter particle, and to our opinion this is
because most dark matter searches focused on mass scales for the
weakly interacting massive particle (WIMP) of the order up to GeV
or even hundred GeV \cite{Oikonomou:2006mh}. However, the
scientific community searching for WIMPs now focuses on mass
scales of eV or much more smaller mass scales. One of the most
promising class of WIMPs with tiny mass is the axion, or any axion
like particle predicted from low scale compactifications of string
theory
\cite{Marsh:2015xka,Marsh:2017yvc,Odintsov:2019mlf,Cicoli:2019ulk,Fukunaga:2019unq,Caputo:2019joi,maxim,Auriol:2018ovo,Ioannisian:2017srr}.
In fact several experimental and observational proposals exist
already in the literature
\cite{Du:2018uak,Henning:2018ogd,Ouellet:2018beu,Safdi:2018oeu,Rozner:2019gba,Avignone:2018zpw,Caputo:2018vmy,Caputo:2018ljp,Millar:2016cjp,Majorovits:2017ppy,Anastassopoulos:2017ftl,TheMADMAXWorkingGroup:2016hpc},
see also \cite{Lawson:2019brd}, and many researchers aim to
identify in the laboratory this elusive particle. To our opinion,
if the axion exists, it will be found in the next 10-15 years due
to extensive experimental searches now performed, mainly based on
the fact that axions and photons interact in the presence of
magnetic fields
\cite{Balakin:2009rg,Balakin:2012up,Balakin:2014oya}, and the
axion is the last resort of particle dark matter, unless nature
hides the WIMPs to a supersymmetry breaking scale, and this will
be another surprise for theorists. The axion is known to provide a
very good candidate for dark matter, or at least some of the dark
matter existing in the Universe \cite{Marsh:2015xka}, and when we
refer to the axion, this should not be confused with the QCD
axion, which is a pseudo-scalar field, but we refer to any
axion-like particle, with a primordial $U(1)$ broken symmetry.
Particularly, these are known as misalignment axion fields, which
are canonical scalar fields. In Ref.~\cite{Odintsov:2019evb}, we
investigated an effective model of $F(R)$ gravity in the presence
of an misalignment axion field, and by using the axion field
dynamics known from the literature \cite{Marsh:2015xka}, we
demonstrated that the axion via a non-minimal coupling to the
scalar curvature affects the late-time era, making possible for
the dark energy era to be realized. Some striking features of the
model were the facts that firstly, the early-time era was
controlled by the well-known $R^2$ gravity
\cite{Starobinsky:1982ee}, and the axion field after it starts
oscillating, for cosmic times for which $m_a\geq H$, so during and
after the reheating era, its energy density scales as $\rho_a\sim
a^{-3}$. Effectively it describes a dark matter condensate field,
which on average yields an effective equation of state (EoS)
parameter $\langle w \rangle =0$.

In this work, we shall further explore several theoretical and
phenomenological features of the model \cite{Odintsov:2019evb},
and of some extensions of this model, to include non-canonical
scalar fields. So the study is focused on theoretical problems and
implications of $F(R)$ gravity in the presence of a non-canonical
in general scalar field, and we emphasize on the phenomenology of
the theory containing the canonical scalar field, which is a
subcase. Particularly, we shall consider the Einstein frame
picture of the scalar-$F(R)$ gravity, and we shall derive the
no-ghost conditions. Also we shall consider the general dynamics
of the model and we shall demonstrate that the resulting
gravitational equations of motion can be used as a reconstruction
technique, which enables us to realize an arbitrary cosmological
evolution. In addition, we shall study the gravitational wave
spectrum, and we shall find several constraints on the
scalar-$F(R)$ gravity theory. Finally, the energy conditions are
considered, and by using these, for the radiation domination era,
we evince how the scalar-$F(R)$ gravity theory can be constrained,
and specifically how the non-minimal coupling of the axion to the
curvature can be constrained.

The motivation for using extensions of $F(R)$ gravity including
scalar fields is threefold. Firstly, $F(R)$ gravity has been
considered and actively investigated as a model quite elegantly
describing dark energy \cite{Capozziello:2003gx,Nojiri:2003ft}. In
the $F(R)$ gravity, appears a scalar particle called scalaron
which interacts with the standard model particles very weakly and
therefore the lifetime of this particle could be rather long. Then
some scenarios where the scalaron can be dark matter particle have
been also considered
\cite{Nojiri:2008nt,Cembranos:2008gj,Katsuragawa:2016yir}. This
scenario is very fascinating because the $F(R)$ gravity might
explain both of the dark energy and dark matter. Since the
scalaron is directly interacting with the standard model
particles, however, the scalaron can decay into the standard
particles even if the interaction is very weak and we need to
consider some scenario to suppress this decay. One easy way to
avoid this problem is to introduce another scalar field and we may
regard the particle corresponding to the field as a dark matter.
Secondly, the Starobinsky model  seems to be quite robust against
the Planck 2018 \cite{Akrami:2018odb} constraints on inflation,
and thirdly, after the Higgs particle discovery
\cite{Aad:2012tfa}, the scalar fields seem to be an inherent
constituent of the primordial era. Thus in this work we shall
further explore the theoretical and phenomenological implications
of the misalignment axion-$F(R)$ gravity.

This paper is organized as follows: In section II, we briefly
recall the essential features of the misalignment axion-$F(R)$
gravity model developed in Ref.~\cite{Odintsov:2019evb}. In
section III, we study the Einstein frame picture of the model, and
we derive the no-ghost criteria, while in section IV, we discuss
the dynamics of the model. We also provide a general
reconstruction technique that can be used to realize any arbitrary
cosmological evolution. In section V we study the gravitational
waves of the axion-$F(R)$ gravity model, and we consider several
constraints related with the observations. In section VI, we
discuss the energy conditions, and we study the dynamics of the
model for several eras, and we demonstrate how the radiation
domination era may impose a constraint on the non-minimal coupling
of the axion to the scalar curvature. Finally, the conclusions
follow at the end of the paper.

\section{Overview of the non-minimally Coupled Axion-$F(R)$ Gravity Model\label{seca}}

In this section we shall briefly recall the main outcomes of the
model developed in Ref.~\cite{Odintsov:2019evb}, in order to put
the present paper in proper context, and to have a reference point
for the studies that follow. The axion-like $F(R)$ gravity model
of Ref.~\cite{Odintsov:2019evb} was based on the following
gravitational action,
\begin{equation}
\label{mainaction} \mathcal{S}=\int d^4x\sqrt{-g}\left[
\frac{1}{2\kappa^2}F(R)+\frac{1}{2\kappa^2}h(\phi)G(R)-\frac{1}{2}\partial^{\mu}\phi\partial_{\mu}\phi-V(\phi)
\right]\, ,
\end{equation}
with $\kappa^2=\frac{1}{8\pi G}$, and $G$ stands for Newton's
gravitational constant. The $F(R)$ gravity chosen in Ref.~\cite{Odintsov:2019evb}
was the well-known $R^2$ gravity,
\begin{equation}\label{starobinsky}
F(R)=R+\frac{1}{36H_i^2}R^2\, .
\end{equation}
Also the coupling $h(\phi)$ to $G(R)$, was assumed to be,
\begin{equation}\label{hphichoice}
h(\phi)\sim\frac{1}{\phi^{\delta}}\, ,
\end{equation}
choosing $\delta>0$, and also the $G(R)$ function is,
\begin{equation}\label{GRfunction}
G(R)\sim R^{\gamma}\, ,
\end{equation}
with the parameter $\gamma$ belonging to the interval $0<\gamma
<0.75$. The model of Ref.~\cite{Odintsov:2019evb} is an effective
model which relies mostly on the axion field dynamics. The axion
field at early times, when the axion mass $m_a$ satisfies $m_a\ll
H$, is considered frozen in its vacuum expectation value $\phi_i$,
with the potential being,
\begin{equation}\label{axionpotential}
V(\phi(t))\simeq \frac{1}{2}m_a^2\phi^2_i(t)\, .
\end{equation}
The following initial conditions at early times are assumed,
\begin{equation}\label{axioninitialconditions}
\dot{\phi}(t_i)=\zeta \ll 1,\,\,\,\phi(t_i)=f_a\theta_a\, ,
\end{equation}
with $t_i$ being the cosmic time corresponding to the inflationary
era, and $f_a$ being the axion decay constant. In view of the
initial conditions, the axion is considered frozen, thus it
contributes a minor overdamped cosmological constant during
inflation, therefore the evolution during the inflationary era is
governed solely by the $R^2$ gravity. As the Universe expands and
cools, when $m_a\geq H$, the axion field starts to oscillate. For
all cosmic times, the axion field equation of motion is,
\begin{equation}\label{scalarfieldeqnduringradiation}
\ddot{\phi}+3H\dot{\phi}+m_a^2\phi=0\, ,
\end{equation}
so for $m_a\geq H$, by assuming a slow-varying oscillatory
behavior for the axion, of the form,
\begin{equation}\label{solutionaxionradandaft}
\phi (t)= A(t)\cos (m_a t)\, ,
\end{equation}
with $A(t)$ quantifying the slow-varying behavior, since it
satisfies,
\begin{equation}\label{Atsolutionconstraints}
\frac{\dot{A}}{m_a}\sim \frac{H}{m_a}\sim \epsilon\ll 1\, ,
\end{equation}
by solving the equation of motion~(\ref{scalarfieldeqnduringradiation})
with the ansatz~(\ref{solutionaxionradandaft}), we get the solution,
\begin{equation}\label{solutionforA}
A\sim a^{-3/2}\, .
\end{equation}
Hence, for cosmic times after the inflationary era, during the
reheating and until late times, the scalar field behaves as,
\begin{equation}\label{solutionaxionradandaftpaper2}
\phi (t)= a^{-3/2}\cos (m_a t)\, .
\end{equation}
Thus, at early times, the axion field has an effective equation of
state (EoS) parameter $w_a=-1$, while for cosmic times $m_a\gg H$,
it has an average EoS of the form $w_a=0$, which describes dark
matter. Finally, as it was shown in Ref.~\cite{Odintsov:2019evb},
the late-time dynamics of the model is governed solely by the term
$\sim h(\phi) G(R)$, with $h(\phi)\sim \phi^{-\delta}$ and
$G(R)\sim R^{\gamma}$ and $0<\gamma <0.75$.

In the following sections we shall discuss several issues related
with the phenomenological aspects of the axion-$F(R)$ gravity
model. Particularly, we shall consider the Einstein frame
counterpart theory, and we discuss the constraints imposed by the
non-ghost condition, and several other theoretical issues related
to the Einstein frame picture. Also we shall present the
reconstruction techniques for the general scalar-$F(R)$ gravity,
which enables us to realize any cosmological solution, given the
Hubble rate. This is particularly useful for cosmological eras in
between the inflationary and dark energy era. Also we shall
discuss the gravitational waves solutions corresponding to early
times, and finally, we shall examine if the energy conditions are
satisfied, and how the radiation domination era may constrain the
coupling of the scalar field $h(\phi)$ to the $G(R)$ gravity term.

\section{The Einstein Frame Picture: Constraints and Limitations \label{Sec1}}

Let us firstly consider the Einstein frame picture of the
scalar-$F(R)$ gravity theory, and we discuss in some detail some
theoretical implications and constraints of the model.

Let us consider a more general gravitational action in comparison
to the one appearing in Eq.~(\ref{mainaction}), of the form,
\begin{equation}
\label{FhGV1} S = \int d^4 x \sqrt{-g} \left[ \frac{1}{2\kappa^2}
\left\{ F(R) + h(\phi) G(R)\right\}
 - \frac{1}{2} \omega (\phi ) \partial_\mu \phi \partial^\mu \phi - V (\phi)
+ \mathcal{L}_\mathrm{matter} \left( g_{\mu\nu}, \Phi_i \right)
\right] \, ,
\end{equation}
where $F(R)$ and $G(R)$ are functions of the scalar curvature $R$
and $h(\phi)$, $\omega(\phi)$, and $V(\phi)$ are functions of the
scalar field $\phi$. Also, $\mathcal{L}_\mathrm{matter} \left(
g_{\mu\nu}, \Phi_i \right) $ is the Lagrangian density of the
matter fluids present -if any- and $\Phi_i$'s express the matter
fields. We now rewrite the action~(\ref{FhGV1}) by introducing the
auxiliary fields $A$ and $B$ as follows,
\begin{equation}
\label{FhGV2}
S_{AB} = \int d^4 x \sqrt{-g} \left[ \frac{1}{2\kappa^2} \left\{ B \left( R - A \right)
+ F(A) + h(\phi) G(A)
 - \frac{1}{2} \omega (\phi ) \partial_\mu \phi \partial^\mu \phi - V (\phi) \right\}
+ \mathcal{L}_\mathrm{matter} \left( g_{\mu\nu}, \Phi_i \right) \right] \, .
\end{equation}
Then by varying the action with respect to the auxiliary scalar
$A$, we obtain,
\begin{equation}
\label{FhGV3}
B = F'(A) + h(\phi) G'(A)\, .
\end{equation}
Then the condition for the absence of the anti-gravity or the
condition that the graviton is not ghost is given by,
\begin{equation}
\label{FhGV4}
F'(A) + h(\phi) G'(A) = F'(R) + h(\phi) G'(R) > 0 \, .
\end{equation}
Obviously, this condition is satisfied by the gravitational action
(\ref{mainaction}). Upon redefining the scalar field $B$ by using
a new scalar field $\sigma$ as $B=\e^{\sigma}$, we assume that
Eq.~(\ref{FhGV3}) can be solved with respect to $A$ as $A=A(\phi,
\sigma)$. Then the action~(\ref{FhGV2}) can be rewritten as
follows,
\begin{align}
\label{FhGV5}
S_{\sigma\phi} = \int d^4 x \sqrt{-g} & \left[ \frac{1}{2\kappa^2} \left\{ \e^\sigma
\left( R - A\left( \phi , \sigma\right) \right)
+ F\left(A \left( \phi , \sigma\right) \right)
+ h(\phi) G \left( A\left( \phi , \sigma\right) \right)
 - \frac{1}{2} \omega (\phi ) \partial_\mu \phi \partial^\mu \phi - V (\phi) \right\}
\right. \nn
& + \mathcal{L}_\mathrm{matter} \left( g_{\mu\nu}, \Phi_i \right) \Bigr] \, .
\end{align}
By the scale transformation of the metric,
\begin{equation}
\label{FhGV6}
g_{\mu\nu} = \e^{-\sigma} {\tilde g}_{\mu\nu} \, ,
\end{equation}
the action~(\ref{FhGV5}) can be rewritten in the Einstein frame,
and it is equal to,
\begin{align}
\label{FhGV7}
S_\mathrm{E} =& \int d^4 x \sqrt{-\tilde g} \left[ \frac{1}{2\kappa^2} \left\{
\tilde R - \frac{3}{2} \partial_\mu \sigma \partial^\mu \sigma
 - \frac{1}{2} \e^{-\sigma} \omega (\phi ) \partial_\mu \phi \partial^\mu \phi
 - U \left( \phi, \sigma \right) \right\}
+ \e^{-2\sigma} \mathcal{L}_\mathrm{matter}
\left( \e^{-\sigma} {\tilde g}_{\mu\nu}, \Phi_i \right) \right] \, , \nn
U \left( \phi, \sigma \right) \equiv &
\e^{-\sigma} A\left( \phi , \sigma\right)
+ \e^{-2\sigma} \left\{ - F\left(A \left( \phi , \sigma\right) \right)
 - h(\phi) G \left( A\left( \phi , \sigma\right) \right)
+ V (\phi) \right\} \, .
\end{align}
In addition to Eq.~(\ref{FhGV4}), the condition that ensures the
absence of ghost degrees of freedom is,
\begin{equation}
\label{FhGV8} \omega (\phi ) > 0 \, ,
\end{equation}
which is clearly satisfied by the gravitational action~(\ref{mainaction}).
Also, let us quote here a marginal remark,
related to the swampland criteria, which is the following, the
constraint given by the so-called swampland conjecture could have
the following form,
\begin{equation}
\label{FhGV9}
\sqrt{ \frac{1}{3} \left( \frac{\partial U}{\partial \sigma} \right)^2
+ \frac{\e^\sigma}{\omega(\phi)} \left( \frac{\partial U}{\partial \phi} \right)^2 }
\geq c U \, ,
\end{equation}
with a constant $c$. Since $B=\e^\sigma$, Eq.~(\ref{FhGV3})
yields,
\begin{equation}
\label{FhGV9B}
0=- \e^\sigma d\sigma + \left( F''(A) + h(\phi) G''(A) \right) dA + h'(\phi) G'(A) d\phi\, ,
\end{equation}
and we obtain,
\begin{equation}
\label{FhGV9C}
\frac{\partial A(\phi,\sigma)}{\partial \phi}
= - \frac{h'(\phi) G'(A)}{F''(A) + h(\phi) G''(A)}\, , \quad
\frac{\partial A(\phi,\sigma)}{\partial \sigma}
= - \frac{\e^\sigma}{F''(A) + h(\phi) G''(A)}\, ,
\end{equation}
which in conjunction with Eqs.~(\ref{FhGV3}) and (\ref{FhGV7}),
indicates that,
\begin{align}
\label{FhGV9D}
\frac{\partial U}{\partial \sigma} =
 - \e^{-\sigma} A - 2 \e^{-2\sigma} \left( - F(A) - h(\phi) G(A) + V (\phi) \right) \, , \quad
\frac{\partial U}{\partial \phi} = \e^{-2\sigma} V'(\phi) \, .
\end{align}
Then since $A=R$ in the original Jordan frame action, and
$B=\e^\sigma$, by using Eq.~(\ref{FhGV3}), we can rewrite the
condition (\ref{FhGV9}), as follows,
\begin{align}
\label{FhGV9F}
& \left[ \frac{1}{3}
\left\{ - \frac{R}{F'(R) + h(\phi) G'(R)}
 - \frac{ 2 \left( - F(R) - h(\phi) G(R) + V (\phi) \right)}
{ \left( F'(R) + h(\phi) G'(R) \right)^2} \right\}^2
+ \frac{\omega(\phi) \left( V' (\phi) \right)^2}{\left( F'(R) + h(\phi) G'(R) \right)^3}
\right]^{\frac{1}{2}} \nn
& \geq c \left[ \frac{R}{F'(R) + h(\phi) G'(R)}
+ \frac{ - F(R) - h(\phi) G(R) + V (\phi) }
{ \left( F'(R) + h(\phi) G'(R) \right)^2} \right] \, .
\end{align}
As an example, we consider the following model,
\begin{equation}
\label{Ex1}
F(R) = R + f_0 R^2 \, , \quad h(\phi) G(R) = h_0 \phi^{-\delta} R^\gamma \, .
\end{equation}
Here $f_0$, $h_0$, $\delta$, and $\gamma$ are constants and we may assume
$\delta>0$ and $0<\gamma<3/4$.
Then Eq.~(\ref{FhGV3}) has the following form
\begin{equation}
\label{Ex2}
B = 1 + 2 f_0 R + \gamma h_0 \phi^{-\delta} R^{\gamma - 1}\, ,
\end{equation}
and the condition~(\ref{FhGV9F}) gives
\begin{align}
\label{Ex2} & \sqrt{ \frac{1}{3} \left( R + \left( 2 - \gamma
\right) h_0 \phi^{-\delta} R^\gamma - V(\phi) \right)^2 +
\omega(\phi) \left( V'(\phi) \right)^2 \left( 1 + 2 f_0 R + \gamma
h_0 \phi^{-\delta} R^{\gamma - 1} \right) } \nn & \leq c \left(
f_0 R^2 + \left( \gamma - 1 \right) h_0 \phi^{-\delta} R^\gamma +
V(\phi) \right) \, ,
\end{align}
which constrains the value of $h_0$ and the non-minimal coupling
$h(\phi)$ in general.

\section{General Dynamics in the Jordan Frame: Inflation and Intermediate Eras\label{Sec2}}

In this section we shall discuss the dynamical evolution of the
scalar-$F(R)$ gravity gravitational system, which is governed by
the equations of motion. Also we shall show how the gravitational
equations can be used as a reconstruction technique, and in effect
an arbitrary cosmological evolution may be realized by the theory
at hand.

By varying the action~(\ref{FhGV1}) with respect to the metric, we
obtain the following equations of motion,
\begin{align}
\label{FhGV10}
0 = & \frac{1}{2}g_{\mu\nu} \left( F(R) + h(\phi) G(R) \right)
 - R_{\mu\nu} \left( F'(R) + h(\phi) G'(R) \right)
 - g_{\mu\nu} \Box \left( F'(R) + h(\phi) G'(R) \right)
+ \nabla_\mu \nabla_\nu \left( F'(R) + h(\phi) G'(R) \right) \nn
& + \frac{1}{2} \omega (\phi ) \partial_\mu \phi \partial_\nu \phi
+ \frac{1}{2} g_{\mu\nu} \left( - \frac{1}{2} \omega (\phi ) \partial_\rho \phi \partial^\rho \phi
 - V \left( \phi \right) \right)
+ \frac{\kappa^2}{2}T_{\mathrm{matter}\, \mu\nu}\, ,
\end{align}
where $T_{\mathrm{matter}\, \mu\nu}$ is the energy-momentum tensor
of the matter fluids present. On the other hand, the variation of
the action~(\ref{FhGV1}) with respect to the scalar field $\phi$
gives,
\begin{equation}
\label{FhGV11}
0 = \partial^\mu \left( \omega(\phi) \partial_\mu \phi \right)
+ h'(\phi) G(R) - V'(\phi) \, .
\end{equation}
In a spatially flat FRW universe,
\begin{equation}
\label{JGRG14}
ds^2 = - dt^2 + a(t)^2 \sum_{i=1,2,3} \left(dx^i\right)^2\, ,
\end{equation}
and by assuming that the scalar field $\phi$ depends only on the
cosmological time $t$, the $(t,t)$ and $(i,j)$ components of
(\ref{FhGV10}) have the following form,
\begin{align}
\label{FhGV12}
0 =& -\frac{1}{2}\left( F(R) + h(\phi) G(R) \right)
+ 3\left(H^2 + \dot H\right) \left( F'(R) + h(\phi) G'(R) \right)
 - 3H \frac{d}{dt} \left( F'(R) + h(\phi) G'(R) \right) \nn
& + \frac{1}{4} \omega (\phi ) {\dot\phi}^2 + \frac{1}{2} V \left( \phi \right)
+ \kappa^2 \rho \, ,\\
\label{FhGV13}
0 =& \frac{1}{2} \left( F(R) + h(\phi) G(R) \right)
 - \left(\dot H + 3H^2\right) \left( F'(R) + h(\phi) G'(R) \right)
+ \left( \frac{d^2}{dt^2} + 2 H \frac{d}{dt} \right)
\left( F'(R) + h(\phi) G'(R) \right)
 \nn
& + \frac{1}{4} \omega (\phi ) {\dot\phi}^2 - \frac{1}{2} V \left(
\phi \right) + \kappa^2 p \, ,
\end{align}
where, the Hubble rate $H$ is defined by $H=\dot a/a$ and the
scalar curvature $R$ is given by $R=12H^2 + 6\dot H$. Furthermore
$\rho$ and $p$ are the energy density and the pressure of the
matter fluids, respectively. By using the ambiguity for the
redefinition of the scalar field $\phi$, we may identify $\phi=t$.
Then Eqs.~(\ref{FhGV12}) and (\ref{FhGV13}) can be rewritten as,
\begin{align}
\label{FhGV14}
\omega (\phi ) =& - 2 \left. \left\{ 2 \dot H \left( F'(R) + h(\phi) G'(R) \right)
+ \left( \frac{d^2}{dt^2} - H \frac{d}{dt} \right)
\left( F'(R) + h(\phi) G'(R) \right)
+ \kappa^2 \left( \rho + p \right) \right\} \right|_{t=\phi} \, ,\\
\label{FhGV15}
V \left( \phi \right) =& \Bigl[ F(R) + h(\phi) G(R)
 - \left(4 \dot H + 6 H^2\right) \left( F'(R) + h(\phi) G'(R) \right)
+ \left( \frac{d^2}{dt^2} + 5 H \frac{d}{dt} \right)
\left( F'(R) + h(\phi) G'(R) \right) \nn
& \left. - \kappa^2 \left( \rho - p \right) \Bigr] \right|_{t=\phi}\, ,
\end{align}
Then given an arbitrary cosmological evolution $H(t)$, by using
the equation of state and the continuity equation, we may find the
explicit $t$ dependencies of $\rho $ and $p $. Then the right hand
side of both Eqs.~(\ref{FhGV14}) and (\ref{FhGV15}) is expressed
as a function of the cosmic time $t$. Since $\phi=t$, by replacing
$t$ with $\phi$ in the right hand side of both Eqs.~(\ref{FhGV14})
and (\ref{FhGV15}), we obtain the explicit forms of $\omega(\phi)$
and $V(\phi)$ as functions of the scalar field $\phi$. This
indicates that an arbitrary expansion history of the Universe
described by $H=H(t)$ can be realized by choosing $\omega(\phi)$
and $V(\phi)$ to satisfy Eqs.~(\ref{FhGV14}) and (\ref{FhGV15})
for any form of $F(R)$, $G(R)$, and $h(\phi)$. As an example, we
consider the model in which the scale factor is,
\begin{equation}
\label{FhGVex1} a = \mathcal{A} (t) \equiv a_0 \e^{ \frac{\left(
H_I - H_L \right) t}{1 + H_T t} + H_L t }\, ,
\end{equation}
where $a_0$, $H_I$, $H_T$, and $H_L$ are dimensionful constants.
Then, Eq.~(\ref{FhGVex1}) yields,
\begin{equation}
\label{FhGVex2}
H = \mathcal{H}(t) \equiv \frac{\left( H_I - H_L \right) }{ \left(1 + H_T t\right)^2 } + H_L \, .
\end{equation}
In effect, cosmic times for which $t\ll \frac{1}{H_T}$, $H\sim
H_I$, correspond to the inflationary era in the early Universe and
cosmic times for which $t\gg \frac{1}{H_T}$, $H\sim H_L$,
correspond to the accelerated expansion of the late Universe. We
also find that the scalar curvature $R$ is given by,
\begin{equation}
\label{FhGVex3}
R = \mathcal{R} (t)
\equiv 12 \left\{ \left( \frac{\left( H_I - H_L \right) }{ \left(1 + H_T t\right)^2 }
+ H_L \right)^2 - \frac{\left( H_I - H_L \right) H_T }{ \left(1 + H_T t\right)^3 }
\right\} \, .
\end{equation}
When the matter has a constant EoS parameter $w\equiv
\frac{p}{\rho}$, $\rho$ is given by $\rho = \rho_0 a^{-3 \left(
1+w \right)}$ with a constant $\rho_0$. Then Eqs.~(\ref{FhGV14})
and (\ref{FhGV15}) indicate that the Universe whose evolution in
the expansion is given by (\ref{FhGVex1}) is realized by choosing,
\begin{align}
\label{FhGV14ex}
\omega (\phi ) =& - 4 \mathcal{H}' \left(\phi \right)
\left( F'\left(\mathcal{R} \left(\phi \right) \right)
+ h(\phi) G'\left(\mathcal{R}\left(\phi\right) \right) \right)
 - 2\left( \frac{d^2}{d\phi^2} - \mathcal{H} \left( \phi \right) \frac{d}{d\phi} \right)
\left( F'\left( \mathcal{R} \left( \phi \right) \right)
+ h(\phi) G' \left( \mathcal{R} \left( \phi \right) \right) \right) \nn
& -2 \kappa^2 \left( 1 + w \right) \rho_0
\left( \mathcal{A} \left(\phi \right) \right)^{-3(1+w)} \, ,\\
\label{FhGV15ex}
V \left( \phi \right) =& F \left( \mathcal{R} \left( \phi \right) \right)
+ h(\phi) G \left( \mathcal{R} \left( \phi \right) \right)
 - \left(4 \mathcal{H}' \left( \phi \right) + 6 \mathcal{H} \left(\phi\right)^2\right)
\left( F' \left( \mathcal{R} \left( \phi \right) \right)
+ h(\phi) G' \left( \mathcal{R} \left( \phi \right) \right) \right) \nn
& + \left( \frac{d^2}{d\phi^2} + 5 \mathcal{H} \left( \phi \right) \frac{d}{d\phi} \right)
\left( F' \left( \mathcal{R} \left( \phi \right) \right)
+ h(\phi) G' \left( \mathcal{R} \left( \phi \right) \right) \right)
 - \kappa^2 \left( 1 - w \right) \rho_0 \left( \mathcal{A} \left(\phi \right)
\right)^{-3(1+w)}  \, .
\end{align}
Here the functions $F(R)$, $G(R)$, and $h(\phi)$ can be arbitrary but we may choose
them as in (\ref{Ex1}).

\section{Gravitational Waves of the Axion-$F(R)$ Gravity Theory\label{Sec3}}

We now consider the gravitational wave based on the action~(\ref{FhGV7})
in the Einstein frame by considering the
perturbation of the background ${\tilde g}_{\mu\nu} = {\tilde
g^{(0)}}_{\mu\nu}$ as ${\tilde g}_{\mu\nu} = {\tilde
g^{(0)}}_{\mu\nu} + {\tilde h}_{\mu\nu}$ in the Einstein equation,
\begin{equation}
\label{FhGV16}
{\tilde R}_{\mu\nu} - \frac{1}{2} {\tilde g}_{\mu\nu} {\tilde R}
= 3 \partial_\mu \sigma \partial_\nu \sigma
+ \e^{-\sigma} \omega (\phi ) \partial_\mu \phi \partial_\nu \phi
+ {\tilde g}_{\mu\nu} \left( - \frac{3}{2} \partial_\rho \sigma \partial^\rho \sigma
 - \frac{1}{2} \e^{-\sigma} \omega (\phi ) \partial_\rho \phi \partial^\rho \phi
 - U \left( \phi, \sigma \right) \right)
+ \kappa^2 {\tilde T}_{\mathrm{matter}\, \mu\nu} \, .
\end{equation}
Here the matter energy momentum tensor ${\tilde
T}_{\mathrm{matter}}^{\mu\nu}$ in the Einstein frame is,
\begin{equation}
\label{FhGV17O}
{\tilde T}_{\mathrm{matter}}^{\mu\nu} \equiv \frac{2}{\sqrt{- \tilde g}}
\frac{\partial \left( \sqrt{- \tilde g} \e^{-2\sigma} \mathcal{L}_\mathrm{matter}
\left( \e^{-\sigma} {\tilde g}_{\mu\nu}, \Phi_i \right) \right)}
{\partial {\tilde g}_{\mu\nu}} \, .
\end{equation}
If the matter fluids couple minimally with gravity, that is, if
the matter Lagrangian $ \mathcal{L}_\mathrm{matter} \left(
\e^{-\sigma} {\tilde g}_{\mu\nu}, \Phi_i \right) =
\mathcal{L}_\mathrm{matter} \left( g_{\mu\nu}, \Phi_i \right)$
does not include any derivative of the metric $g_{\mu\nu}$, the
matter energy momentum tensor ${\tilde
T}_{\mathrm{matter}}^{\mu\nu}$ in the Einstein frame is related
with the matter energy momentum tensor
$T_{\mathrm{matter}}^{\mu\nu}$ in the original Jordan frame as
${\tilde T}_{\mathrm{matter}}^{\mu\nu} = \e^{-3\sigma}
T_{\mathrm{matter}}^{\mu\nu}$, that is, ${\tilde
T}_{\mathrm{matter}\, \mu\nu} = \e^{-\sigma} T_{\mathrm{matter}\,
\mu\nu}$. When we consider the gravitational wave, we often use
the transverse and traceless gauge conditions. Since we are
considering the scale transformation~(\ref{FhGV6}), if
$h_{\mu\nu}$, which is defined by the fluctuation from the
background metric $g_{\mu\nu} = g^{(0)}_{\mu\nu}$ as $g_{\mu\nu} =
g^{(0)}_{\mu\nu} + h_{\mu\nu}$ in the original frame in
the action~(\ref{FhGV1}), satisfies the transverse and traceless gauge
conditions,
\begin{equation}
\label{FhGV17} \nabla^\mu h_{\mu\nu} = g^{(0)\, \mu\nu} h_{\mu\nu}
= 0 \, .
\end{equation}
However, the scale transformed fluctuation ${\tilde h}_{\mu\nu} =
\e^\sigma h_{\mu\nu}$ does not always satisfies the first
condition in (\ref{FhGV17}), although the second condition is
trivially satisfied, ${\tilde g}^{(0)\, \mu\nu} {\tilde
h}_{\mu\nu} = \e^{-\sigma} g^{(0)\, \mu\nu} \e^\sigma h_{\mu\nu} =
g^{(0)\, \mu\nu} h_{\mu\nu} = 0$. For the first condition in
(\ref{FhGV17}), under the scale transformation, we find,
\begin{align}
\label{FhGV18}
& \tilde{\nabla}^\mu \tilde{h}_{\mu}^{\ \nu}
= \e^{-\sigma} \nabla^\mu h_{\mu\nu}
+ 4 \e^{-\sigma} g^{(0)\, \mu\tau} g^{(0)\, \nu\rho} \sigma_{,\tau} h_{\mu\rho}
 - \e^{-\sigma} g^{(0)\, \nu\rho}\sigma_{,\rho} g^{(0)\, \mu\tau} h_{\mu\tau}
= 4 \e^{-\sigma} g^{(0)\, \mu\tau} g^{(0)\, \nu\rho} \sigma_{,\tau} h_{\mu\rho} \, .
\end{align}
Then if we assume that a homogeneous and isometric background
metric, and therefore $\sigma$ only depends on the cosmological
time $t$ and also $g^{(0)}_{ti}=0$, then if we consider the
perturbation with $h_{t\mu}=0$ since we are considering the
massless spin 2 mode, we find,
\begin{equation}
\label{FhGV19} \tilde{\nabla}^\mu \tilde{h}_{\mu}^{\ \nu} =
{\tilde g}^{(0)\, \mu\nu} {\tilde h}_{\mu\nu} = 0\, .
\end{equation}
Therefore the gauge conditions in (\ref{FhGV17}) for the graviton are not changed
by the scale transformation,
\begin{equation}
\label{FhGV21}
{\tilde\nabla}^\mu {\tilde h}_{\mu\nu} = {\tilde g}^{(0)\, \mu\nu} {\tilde h}_{\mu\nu} = 0 \, .
\end{equation}
Then under the condition~(\ref{FhGV21}), the equation for the
gravitational wave can be written as follows,
\begin{align}
\label{FhGV22}
0 =& \frac{1}{2\kappa^2}\left(- \frac{1}{2}\left( - {\tilde \Box}^{(0)} {\tilde h}_{\mu\nu}
 - 2{\tilde R}^{(0)\, \lambda\ \rho}_{\ \ \ \ \ \nu\ \mu} {\tilde h}_{\lambda\rho}
+ {\tilde R}^{(0)\, \rho}_{\ \ \ \ \ \mu}{\tilde h}_{\rho\nu}
+ {\tilde R}^{(0)\, \rho}_{\ \ \ \ \ \nu}{\tilde h}_{\rho\mu} \right)
+ \frac{1}{2} R^{(0)} h_{\mu\nu}
 - \frac{1}{2}{\tilde g}^{(0)}_{\mu\nu} {\tilde h}_{\rho\sigma}
{\tilde R}^{(0)\, \rho\sigma} \right) \nn
& + {\tilde h}_{\mu\nu} \left( - \frac{3}{2} \partial_\rho \sigma \partial^\rho \sigma
 - \frac{1}{2} \e^{-\sigma} \omega (\phi ) \partial_\rho \phi \partial^\rho \phi
 - U \left( \phi, \sigma \right) \right)
- {\tilde g}^{(0)}_{\mu\nu} \left( \frac{3}{2} \partial^\rho \sigma \partial^\tau \sigma
 - \frac{1}{2} \e^{-\sigma} \omega (\phi ) \partial^\rho \phi \partial^\tau \phi
\right) {\tilde h}_{\rho\tau} \nn
& + \frac{1}{2} \frac{\partial {\tilde T}_{\mathrm{matter}\, \mu\nu}}
{\partial {\tilde g}_{\rho\tau}}
{\tilde h}_{\rho\tau} \, .
\end{align}
We are now interested in the massless spin two mode, which
satisfies,
\begin{equation}
\label{FhGV23}
{\tilde h}_{it}={\tilde h}_{ti}=h_{it}=h_{ti}=0\, , \quad
\sum_{i=1,2,3} {\tilde h}_{ii}=\sum_{i=1,2,3} h_{ii}=0 \, , \quad i=1,2,3\, ,
\quad {\tilde h}_{tt}=h_{tt}=0 \, .
\end{equation}
In the spatially flat FRW universe in the Einstein frame,
\begin{equation}
\label{JGRG14E}
d{\tilde s}^2 \equiv \e^{\sigma} ds^2 = - d{\tilde t}^2
+ {\tilde a} \left( \tilde t \right)^2 \sum_{i=1,2,3} \left(dx^i\right)^2\, ,
\end{equation}
where $d{\tilde t}\equiv \e^{\frac{\sigma}{2}} dt$ and $\tilde a
\left( \tilde t \right) \equiv \e^{\frac{\sigma}{2}} a(t)$, due to
the isometry in the spacial part, we may assume,
\begin{equation}
\label{FhGV24}
\frac{\partial {\tilde T}_{\mathrm{matter}\, tt}}{\partial{\tilde g}_{ij}} \propto \delta^{ij} \, , \quad
\frac{\partial {\tilde T}_{\mathrm{matter}\, tk}}{\partial {\tilde g}_{ij}}
= \frac{\partial {\tilde T}_{\mathrm{matter}\, kt}}{\partial {\tilde g}_{ij}} = 0 \, .
\end{equation}
We may further assume that the matter energy-momentum tensor~(\ref{FhGV17O})
in the Einstein frame has the following form as in the perfect fluid,
\begin{equation}
\label{FhGV25}
{\tilde T}_{\mathrm{matter}\, \mu\nu} =\tilde\rho {\tilde U}_\mu {\tilde U}_\nu
+ \tilde p {\tilde\gamma}_{\mu\nu} \, .
\end{equation}
Here $\left( \tilde U^\mu \right)$ is the four velocity of the matter
fluid and we now assume $\tilde U^0=1$ and $\tilde U^i=0$. In
Eq.~(\ref{FhGV25}), ${\tilde\gamma}_{\mu\nu}$ is the projection
tensor to the spatial directions perpendicular to $\tilde U^\mu$,
\begin{equation}
\label{FhGV26}
{\tilde\gamma}_{\mu\nu}={\tilde g}_{\mu\nu} + {\tilde U}_\mu {\tilde U}_\nu\, .
\end{equation}
We now also assume that the matter fluid minimally couples with
the metric ${\tilde g}_{\mu\nu}$, that is, the coupling between
the matter fluids and the metric does not include the derivative
of the metric. Then, under the perturbation ${\tilde g}_{\mu\nu} =
{\tilde g}^{(0)}_{\mu\nu} + {\tilde h}_{\mu\nu}$, we find,
\begin{equation}
\label{FhGV27}
\delta \tilde \rho = {\tilde \rho}^{\mu\nu} {\tilde h}_{\mu\nu}\, , \quad
\delta \tilde p = {\tilde p}^{\mu\nu} {\tilde h}_{\mu\nu}\, .
\end{equation}
On the other hand, the variation of ${\tilde U}_\mu$ is given by
using the condition ${\tilde U}^\mu {\tilde U}_\mu = -1$, that is,
\begin{equation}
\label{FhGV29}
0 = 2 \left( \delta {\tilde U}^\mu \right)
+ {\tilde U}^\mu {\tilde U}^\nu {\tilde h}_{\mu\nu}
= {\tilde U}^\mu \left( 2 {\tilde g}^{(0)}_{\mu\nu} \delta {\tilde U}^\nu
+ {\tilde h}_{\mu\nu} {\tilde U}^\nu \right)\, ,
\end{equation}
and hence we have,
\begin{equation}
\label{FhGV30}
0 = 2 \left( \delta {\tilde U}^\mu \right)
+ {\tilde U}^\mu {\tilde U}^\nu {\tilde h}_{\mu\nu}
= {\tilde U}^\mu \left( 2 {\tilde g}^{(0)}_{\mu\nu} \delta {\tilde U}^\nu
+ {\tilde h}_{\mu\nu} {\tilde U}^\nu \right)\, ,
\end{equation}
or equivalently,
\begin{equation}
\label{FhGV31} \delta {\tilde U}^\mu = - \frac{1}{2} {\tilde
g}^{(0)\, \mu\rho} \left( {\tilde h}_{\rho\nu} {\tilde U}^\nu +
l_\rho \right)\, .
\end{equation}
Note that the arbitrary vector $l_\mu$ satisfies ${\tilde U}^\mu
l_\mu = 0$, but we choose $l_\mu=0$ by assuming the isometry. Due
to the isometry, we may assume $\rho^{ij}$ and $p^{ij}$ are
proportional to $\delta^{ij}$, Then the $(t,t)$ and $(t,i)$
components of equation for gravitational wave (\ref{FhGV22}) are
trivially satisfied and the $(i,j)$ component is given by,
\begin{align}
\label{FhGV32}
0 =& \frac{1}{2\kappa^2}\left( \frac{1}{2}
\left( - \partial_{\tilde t}^2 {\tilde h}_{ij} + {\tilde a}^{-2} \triangle {\tilde h}_{ij} \right)
+ \left( 3 \frac{d \tilde H}{d\tilde t} + 4 {\tilde H}^2 \right) {\tilde h}_{ij}
+ {\tilde h}_{ij} \left( \frac{3}{2} \left( \frac{d \sigma}{d\tilde t} \right)^2
+ \frac{1}{2} \e^{-\sigma} \omega (\phi )
\left( \frac{d \phi}{d\tilde t} \right)^2
 - U \left( \phi, \sigma \right) \right) \right)
+ \frac{1}{2} \tilde p {\tilde h}_{ij} \nn
=& \frac{1}{2\kappa^2}\left( \frac{1}{2}
\left( - \partial_{\tilde t}^2 {\tilde h}_{ij} + {\tilde a}^{-2} \triangle {\tilde h}_{ij} \right)
+ \left( 3 \frac{d \tilde H}{d\tilde t} + 4 {\tilde H}^2 \right) {\tilde h}_{ij}
+ {\tilde h}_{ij} \left( \frac{3}{2} \left( \frac{d \sigma}{d\tilde t} \right)^2
+ \frac{1}{2} \e^{-2\sigma} \omega (\phi )
 - U \left( \phi, \sigma \right) \right) \right)
+ \frac{1}{2} \tilde p {\tilde h}_{ij} \, .
\end{align}
Here we have used $\left( \frac{d \phi}{d\tilde t} \right)^2
=\left( \frac{dt}{d\tilde t} \right)^2 = \e^{-\sigma}$. We may
rewrite Eq.~(\ref{FhGV32}) in the original Jordan frame. Since,
${\tilde h}_{\mu\nu} = \e^\sigma h_{\mu\nu}$ and $\tilde a =
\e^{\frac{\sigma}{2}} a$, we find,
\begin{align}
\label{FhGV33}
\frac{\partial}{\partial \tilde t} =& \e^{-\frac{\sigma}{2}} \frac{\partial}{\partial t} \, , \quad
\frac{\partial^2}{\partial \tilde t^2} = \e^{-\sigma}
\left( \frac{\partial^2}{\partial t^2 } - \frac{1}{2}\frac{d\sigma}{dt} \frac{\partial}{\partial t}
\right) \, , \nn
\tilde H =& \frac{1}{\tilde a} \frac{d \tilde a}{d\tilde t}
= \e^{-\frac{\sigma}{2}} \left( \frac{1}{2} \frac{d \sigma }{d t} + H \right) \, , \quad
\frac{d \tilde H}{d\tilde t} = \e^{-\tilde \sigma} \left( \frac{1}{2} \frac{d^2 \sigma}{dt^2}
+ \frac{1}{4} \left( \frac{d \sigma}{dt} \right)^2 + \frac{d H}{dt} \right) \, , \nn
\frac{\partial^2 {\tilde h}_{ij}}{\partial \tilde t^2} =&
\frac{\partial^2 h_{ij} }{\partial t^2} + \frac{d^2 \sigma}{dt^2} h_{ij}
+ \left( \frac{d\sigma}{dt} \right)^2 h_{ij}
+ 2 \frac{d \sigma}{dt} \frac{\partial h_{ij}}{\partial t}
 - \frac{1}{2} \left( \frac{d \sigma}{dt} \right)^2 \frac{\partial h_{ij}}{\partial t}
 - \frac{1}{2}\frac{d\sigma}{dt} \frac{\partial h_{ij}}{\partial t} \nn
=& \frac{\partial^2 h_{ij} }{\partial t^2} + \frac{d^2 \sigma}{dt^2} h_{ij}
+ \frac{1}{2} \left( \frac{d\sigma}{dt} \right)^2 h_{ij}
+ \frac{3}{2} \frac{d \sigma}{dt} \frac{\partial h_{ij}}{\partial t} \, ,
\end{align}
Then, due to the fact that $\tilde p = \e^{-2\sigma} p$ from
(\ref{FhGV25}) (${\tilde T}_{\mathrm{matter}\, \mu\nu} =
\e^{-\sigma} T_{\mathrm{matter}\, \mu\nu}$ and
${\tilde\gamma}_{\mu\nu} = \e^\sigma \gamma_{\mu\nu}$),
Eq.~(\ref{FhGV32}) is rewritten as follows,
\begin{align}
\label{FhGV34}
0 =& \frac{1}{2\kappa^2}\left( \frac{1}{2}
\left( - \partial_t^2 h_{ij} + \frac{3}{2} \dot\sigma \partial_t h_{ij}
+ \left( \ddot\sigma + {\dot\sigma}^2 \right) h_{ij} + a^{-2} \triangle h_{ij} \right) \right. \nn
& \left. + \left( 3 \dot H + 4 H^2 + \frac{3}{2} \ddot \sigma + \frac{13}{4} {\dot \sigma}^2
+ 4 \dot\sigma H + \frac{1}{2} \omega (\phi ) {\dot\phi}^2
 - \e^{\sigma} U \left( \phi, \sigma \right) \right) h_{ij} \right)
+ \frac{1}{2}\e^{-\sigma} p h_{ij} \, .
\end{align}
We should note that $\e^\sigma$ is given in Eq.~(\ref{FhGV3}) and
$U \left( \phi, \sigma \right)$ is given in Eq.~(\ref{FhGV7}).
Then in terms of the Jordan frame, we find,
\begin{align}
\label{FhGV35}
\sigma = & \ln \left( F'(R) + h(\phi) G'(R) \right) \, , \nn
U \left( \phi, \sigma \right) = & \e^{-\sigma} R
+ \e^{-2\sigma} \left\{ - F\left( R \right) - h(\phi) G \left( R \right)
+ V (\phi) \right\} \nn
= & \frac{R \left( F'(R) + h(\phi) G'(R) \right)
 - F\left( R \right) - h(\phi) G \left( R \right)
+ V (\phi) }{\left( F'(R) + h(\phi) G'(R) \right)^2} \, .
\end{align}
By choosing $\phi=t$ and using (\ref{FhGV14}) and (\ref{FhGV15}),
we obtain,
\begin{align}
\label{FhGV36}
\frac{1}{2} & \omega (\phi ) {\dot\phi}^2 - \e^{\sigma} U \left( \phi, \sigma \right) \nn
=& - \left( F'(R) + h(\phi) G'(R) \right)^{-1}
\left( 2 \dot H \left( F'(R) + h(\phi) G'(R) \right)
+ \left( \frac{d^2}{dt^2} - H \frac{d}{dt} \right)
\left( F'(R) + h(\phi) G'(R) \right)
+ \kappa^2 \left( \rho + p \right) \right. \nn
& - \left( \left( 12 H^2 + 6 \dot H \right) \left( F'(R) + h(\phi) G'(R) \right)
 - \left(4 \dot H + 6 H^2\right) \left( F'(R) + h(\phi) G'(R) \right) \right. \nn
& \left. \left. + \left( \frac{d^2}{dt^2} + 5 H \frac{d}{dt} \right)
\left( F'(R) + h(\phi) G'(R) \right) - \kappa^2 \left( \rho - p
\right) \right) \right) \nn
=& - \left( F'(R) + h(\phi) G'(R) \right)^{-1}
\left( - 6 \left( H^2 + H \frac{d}{dt} \right) \left( F'(R) + h(\phi) G'(R) \right)
+ 2 \kappa^2 \rho \right) \nn
=& 6 H^2 + 6 H \dot\sigma + 2 \kappa^2 \e^{-\sigma} \rho \, ,
\end{align}
and we can further rewrite Eq.~(\ref{FhGV34}) as follows,
\begin{equation}
\label{FhGV37}
0 = \frac{1}{2\kappa^2}\left( \frac{1}{2}
\left( - \partial_t^2 h_{ij} + \frac{3}{2} \dot\sigma \partial_t h_{ij}
+ a^{-2} \triangle h_{ij} \right)
+ \left( 3 \dot H + 10 H^2 + 2 \ddot \sigma + \frac{17}{4} {\dot \sigma}^2
+ 10 \dot\sigma H \right) h_{ij} \right)
+ \frac{1}{2}\e^{-\sigma} \left( \rho + p \right) h_{ij} \, .
\end{equation}
As an example, we consider the quasi-de Sitter spacetime case,
\begin{equation}
\label{FhGV37}
H=H_0 + H_1 \left( t - t_0 \right)