# Kernel Selection for Gaussian Process in Cosmology: with Approximate
Bayesian Computation Rejection and Nested Sampling

Hao Zhang Institute for Frontiers in Astronomy and Astrophysics, Beijing
Normal University, Beijing 102206, China  Department of Physics, Beijing
Normal University, Beijing 100875, China  Key Laboratory for Particle
Astrophysics, Institute of High Energy Physics, CAS 19B Yuquan Road, Beijing
100049, China  School of Astronomy and Space Science, University of Chinese
Academy of Sciences, 19A Yuquan Road, Beijing 100049, China
<EMAIL_ADDRESS>Yu-Chen Wang Department of Astronomy, School of Physics,
Peking University, Beijing 100871, China  Kavli Institute for Astronomy and
Astrophysics, Peking University, Beijing 100871, China  Tong-Jie Zhang
Institute for Frontiers in Astronomy and Astrophysics, Beijing Normal
University, Beijing 102206, China  Department of Astronomy, Beijing Normal
University, Beijing 100875, China  Institute for Astronomical Science, Dezhou
University, Dezhou 253023, China<EMAIL_ADDRESS>Tingting Zhang  College
of Command and Control Engineering, PLA Army Engineering University, Nanjing
210000, China<EMAIL_ADDRESS>
###### Abstract

Gaussian Process (GP) has gained much attention in cosmology due to its
ability to reconstruct cosmological data in a model-independent manner. In
this study, we compare two methods for GP kernel selection: Approximate
Bayesian Computation (ABC) Rejection and nested sampling. We analyze three
types of data: cosmic Chronometer data (CC), Type Ia Supernovae (SNIa), and
Gamma Ray Burst (GRB), using five kernel functions. To evaluate the
differences between kernel functions, we assess the strength of evidence using
Bayes factors. Our results show that, for ABC Rejection, the Matérn kernel
with $\nu$=5/2 (M52 kernel) outperformes the commonly used Radial Basis
Function (RBF) kernel in approximating all three datasets. Bayes factors
indicate that the M52 kernel typically supports the observed data better than
the RBF kernel, but with no clear advantage over other alternatives. However,
nested sampling gives different results, with the M52 kernel losing its
advantage. Nevertheless, Bayes factors indicate no significant dependence of
the data on each kernel.

Gaussian Processes regression(1930), Astronomy data analysis(1858), Bayes
factor(1919), Observational cosmology(1146), Astrostatistics strategies(1885)

##  1 Introduction

Given the widespread use of Gaussian Process in cosmology, it is crucial to
address various related issues. Cosmologists often use Gaussian Process to
reconstruct cosmological data to anticipate specific cosmological parameters
at a given redshift and constrain these parameters through prediction (Gómez-
Valent & Amendola (2018); Dhawan et al. (2021)). Gaussian Process can also be
utilized to address the Hubble Tension(Mehrabi & Basilakos, 2020), constrain
cosmological parameters in a non-parametric way(Mukherjee & Banerjee, 2021),
and deal with relevant issues in the realm of dark energy (Seikel et al.,
2012), as it can reconstruct data in a model-independent manner. As Gaussian
Process is widely used, the specifics of its calculation, including
hyperparameters, quantity of hyperparameters, and choice of kernels, can
impact the reconstruction of cosmological data and our predictions.

The choice of hyperparameters, the quantity of hyperparameters, and the choice
of kernels can all affect the results of the Gaussian Process. Earlier studies
have addressed many of these issues. For instance, Sun et al. (2021) studied
the influence of the upper bound and lower bound of the two hyperparameters in
the radial basis kernel function (RBF). There were also some articles about
the choice of kernels. However, these studies, such as Shafieloo et al. (2012)
and Hwang et al. (2023), did not provide sufficient details on the choice of
kernel functions and only spent a relatively small portion of their space on
the issue without a clear selection scheme. Others, such as Ó Colgáin &
Sheikh-Jabbari (2021), did not cover a sufficient variety of kernel
functions.Therefore, our research aims to evaluate which kernel outperforms
others for various datasets when various distance functions are used to assess
the effect of reconstruction and the degree to which each kernel is supported
by observed data.

To accomplish our objectives, we use Approximate Bayesian Computation (ABC)
Rejection combined with three distinct distance functions to evaluate the
performance of five kernels against three different cosmological datasets,
including cosmic chronometer data (CC), Type Ia Supernovae (SNIa), and Gamma
Ray Burst data (GRB). Additionally, we calculate the Bayes factor between
every two kernels based on the posterior distribution provided by the ABC
Rejection to show the subtle differences in support strength of our datasets
for the five kernels we test. The Bayes factor $K$ is a powerful tool for
quantitatively comparing the strength of evidence for two models from observed
data given the prior distribution of two different models. This allows us to
determine exactly how superior the ideal kernel function is to alternative
kernels.

This paper’s structure is as follows. In Section 2, we introduce Gaussian
Process (Section 2.1) and other mathematical methods, such as Approximate
Bayesian Computation (ABC) Rejection (Section 2.2), distance function (Section
2.3), Bayes factor (Section 2.4), and nested sampling (Section 2.5), briefly.
In Section 3, we describe the results of calculations using the mathematical
methods mentioned for the cosmological datasets we use, respectively,
involving ABC Rejection, Bayes factor, and nested sampling. In the last
section, Section 4, we summarize our results and discuss the caution that
scientists should exercise in kernel selection when using Gaussian Process to
reconstruct data.

##  2 Methods

###  2.1 Gaussian Process(GP)

The Gaussian Process (GP) is a powerful tool for modeling functions in a
stochastic statistical process. It represents a collection of random variables
with a joint Gaussian distribution from the perspective of function-space
(Williams & Rasmussen, 2006). The observations in GP are presented in a
continuous domain, where a normally distributed random variable is connected
to each point in the input space. By using GP, we can predict the value of a
function at a new point without assuming the form of the function, based on
the observed data points. Specifically, a Gaussian random variable at a
reconstructed point $z$ represents the expected value for GP.

To forecast the function value at $z$, only the covariance function and mean
function of GP are needed. Let $\boldsymbol{Z}=\\{z_{i}\\}$ be a vector of
known redshift data points. The crucial step in obtaining the mean and
covariance function is calculating the covariance matrix K(Z,Z), which is
determined by the observed data, hyperparameters, and kernel function
$k(z_{i},z_{j})$ for arbitrarily chosen reconstructed points $z_{i}$ and
$z_{j}$. The covariance matrix K(Z,Z) displays the relationship between
different reconstructed points, indicating that function values at various
places in GP are not independent of one another. The kernel function returns a
single element, $K_{ij}$, of the covariance matrix K(Z,Z) (Seikel et al.,
2012):

$[\emph{K{(Z,Z)}}]_{ij}=k(z_{i},z_{j}).$ (1)

The primary goal of GP is to determine the mean and variance of the vector
$\boldsymbol{f^{*}}=\\{f(z_{i}^{*})\\}$ containing the reconstructed function
values at prediction points. In order to obtain these values, the mean
$\overline{\boldsymbol{f^{*}}}$ and covariance matrix
cov($\boldsymbol{f^{*}}$) must be calculated using the approach outlined in
Williams & Rasmussen (2006). Specifically, Equation (2) provides the
expressions for these quantities:

$\displaystyle\overline{\boldsymbol{f^{*}}}$
$\displaystyle=K(\boldsymbol{Z^{*}},\boldsymbol{Z})[K(\boldsymbol{Z},\boldsymbol{Z})+C]^{-1}\boldsymbol{y},$
(2) $\displaystyle\mathrm{cov}(\boldsymbol{f^{*}})$
$\displaystyle=K(\boldsymbol{Z^{*}},\boldsymbol{Z^{*}})$
$\displaystyle-K(\boldsymbol{Z^{*}},\boldsymbol{Z})[K(\boldsymbol{Z},\boldsymbol{Z})+C]^{-1}K(\boldsymbol{Z},\boldsymbol{Z^{*}}),$

where $\boldsymbol{Z^{*}}=${$z_{i}^{*}$} represents the vector of redshifts to
be predicted, and $\boldsymbol{y}=\\{y_{i}\\}$ denotes the vector of
observations. The covariance matrix of the observed data is denoted by $C$,
where $C=\text{diag}\\{\sigma_{i}^{2}\\}$ and $\sigma_{i}$ represents the
error value corresponding to each data point. In addition, the calculation of
the covariance matrix depends on the sampling of the hyperparameters
($\sigma_{f}^{2}$,$l$). To prevent the program from encountering a “singular
matrix” error, Cholesky decomposition is used instead of solving the inverse
of the matrix directly throughout the Gaussian process.

The final reconstruction of the GP can be influenced by any variables that
impact the computation of the covariance matrix. This is because the
covariance matrix K(Z,Z) defines both the mean function and covariance
function of the GP. The shape of the kernel function can directly affect the
outcome. In cosmology research, it is commonly assumed that the kernel
function should only depend on the distance between different locations.
Therefore, the RBF kernel is often used due to its simplicity:

  * •

Radial Basis Function (RBF)

$k(z_{i},z_{j})=\sigma_{f}^{2}\exp({-\frac{(z_{i}-z_{j})^{2}}{2l^{2}}}).$ (3)

In Equation (3), the hyperparameters $\sigma_{f}^{2}$ and $l$ are typically
used in kernel functions. However, there are other kernels besides the RBF
kernel that can be used in practical research, such as the Cauchy kernel, the
Matérn($\nu$=5/2) kernel, the Matérn($\nu$=7/2) kernel, and the
Matérn($\nu$=9/2) kernel. These kernels will also be investigated in this
study. Here is a list of the additional kernels:

  * •

Cauchy kernel (CHY)

$\displaystyle
k(z_{i},z_{j})=\sigma_{f}^{2}(1+\frac{(z_{i}-z_{j})^{2}}{2l^{2}})^{-1},$ (4)

  * •

Matérn($\nu$=5/2) kernel (M52)

$\displaystyle k(z_{i},z_{j})$
$\displaystyle=\sigma_{f}^{2}\exp({-\sqrt{5}\frac{|z_{i}-z_{j}|}{l}})$
$\displaystyle(1+\sqrt{5}\frac{|z_{i}-z_{j}|}{l}+5\frac{(z_{i}-z_{j})^{2}}{3l^{2}}),$
(5)

  * •

Matérn($\nu$=7/2) kernel (M72)

$\displaystyle
k(z_{i},z_{j})=\sigma_{f}^{2}\exp({-\sqrt{7}\frac{|z_{i}-z_{j}|}{l}})(1+$
$\displaystyle\sqrt{7}\frac{|z_{i}-z_{j}|}{l}+14\frac{(z_{i}-z_{j})^{2}}{5l^{2}}+7\sqrt{7}\frac{|z_{i}-z_{j}|^{3}}{15l^{3}}),$
(6)

  * •

Matérn($\nu$=9/2) kernel (M92)

$\displaystyle
k(z_{i},z_{j})=\sigma_{f}^{2}\exp({-3\frac{|z_{i}-z_{j}|}{l}})(1+3\frac{|z_{i}-z_{j}|}{l}$
$\displaystyle+27\frac{(z_{i}-z_{j})^{2}}{7l^{2}}+18\frac{|z_{i}-z_{j}|^{3}}{7l^{3}}+27\frac{(z_{i}-z_{j})^{4}}{35l^{4}}).$
(7)

The kernels mentioned earlier are commonly used in various research fields.
However, the M92 kernel is of particular interest as studies have revealed it
to be the most effective kernel for constraining the equation of state
parameter of dark energy (Seikel & Clarkson, 2013), and it is regularly
employed in dark energy surveys. Thus, we have also incorporated the M92
kernel in our study.

The impact of the hyperparameters ($\sigma_{f}^{2}$ and $l$, as explained
earlier) on the computation of the covariance matrix should not be overlooked.
These hyperparameters must be trained for a complete GP in order to determine
their optimal values. In cosmological research, the Log Marginal Likelihood
(LML) is a widely used method for hyperparameter training. The goal of
hyperparameter training is to find the combination of hyperparameters that
correspond to the maximum LML, which is then used as the final set of
hyperparameters to be substituted into the GP for the final result. The
mathematical form of LML is as follows:

$\displaystyle\ln\mathcal{L}$
$\displaystyle=-\frac{1}{2}\boldsymbol{y}^{\top}[K(\boldsymbol{Z},\boldsymbol{Z})+C]^{-1}\boldsymbol{y}$
$\displaystyle-\frac{1}{2}\ln|K(\boldsymbol{Z},\boldsymbol{Z})+C|-\frac{n}{2}\ln
2\pi,$ (8)

where $n$ is the dimension of the vector $\boldsymbol{Z}$. It’s important to
note that other methods can also be used to obtain hyperparameters. When the
LML reaches its maximum value, the relevant hyperparameters yield the most
likely manifestation of the function. In practice, almost all GPs are carried
out by maximizing the LML function.

To illustrate how our data can be used to complete the GP, let’s consider an
example. Suppose we have an $H(z)-z$ dataset with $n$ data points, and we
assume that the hyperparameters ($\sigma_{f}^{2}$ and $l$) are sampled from a
given distribution. With the kernel function chosen, we can construct a kernel
function matrix $K(\boldsymbol{X,X})$ for the $n$ redshift values contained in
the data, where $\boldsymbol{X}$ is a column vector consisting of the
redshifts of the observation points. We also have $n^{\prime}$ prediction
points at which we want the GP to predict the value of the Hubble parameter
$H(z)$, and the vector consisting of these $n^{\prime}$ redshifts is denoted
as $\boldsymbol{X^{*}}$. By combining the error $\sigma_{n}$ at each data
point and the actual observed Hubble parameter $H(z)$ obtained (here we use
$\boldsymbol{y}$ to refer to the column vector made up of observations), we
can start the GP:

$\displaystyle L$ $\displaystyle=\text{cholesky}(K(X,X)+\sigma_{n}^{2}I),$
$\displaystyle\boldsymbol{\alpha}$
$\displaystyle=L^{\top}\backslash(L\backslash\boldsymbol{y}),$ (9)

where $I$ is the identity matrix. It has the same dimension as
$K(\boldsymbol{X,X})$. Now we can get the mean values:

$\overline{\boldsymbol{f^{*}}}=K(\boldsymbol{X^{*},X})\boldsymbol{\alpha}.$
(10)

Next we construct the variance:

$\displaystyle\boldsymbol{v}$ $\displaystyle=L\backslash
K(\boldsymbol{X^{*},X}),$ $\displaystyle\mathrm{cov}(\boldsymbol{f^{*}})$
$\displaystyle=K(\boldsymbol{X^{*},X^{*}})-\boldsymbol{v}^{\top}\boldsymbol{v}.$
(11)

Also, we can get the log-marginal-likelihood:

$\ln\mathcal{L}=-\frac{1}{2}\boldsymbol{y}^{\top}\boldsymbol{\alpha}-\frac{1}{2}\ln|K(\boldsymbol{X},\boldsymbol{X})+\sigma_{n}^{2}I|-\frac{n}{2}\ln
2\pi.$ (12)

It’s important to note that in the previous discussion of the GP principle, we
replace $\sigma_{n}^{2}I$ with $C$. However, this is only a convenient
substitution for presentation purposes. If you need to write an actual GP
algorithm, the steps mentioned above are what you need to follow. This is also
the step where we apply the principles of GP to model the data (Williams &
Rasmussen, 2006).

In this study, we evaluate the effectiveness of various GP kernels using three
datasets that correspond to three different cosmological probes: cosmic
chronometer (CC), type Ia supernovae (SNIa), and Gamma Ray Burst (GRB). The CC
dataset is presented in Table 1, while the SNIa data, obtained from Scolnic et
al. (2018), are binned and consist of 40 data points. The original dataset
contains 1048 points, but the binned version is used in this study. The
covariance matrix, provided by Scolnic et al. (2018), is used directly in the
GP calculations. The GRB dataset, obtained from Demianski et al. (2017), has a
large number of data points in the high-redshift zone, making it an excellent
source of information about the high-redshift universe. We choose not to
transform the GRB and SNIa data into redshift with other parameters since this
process is model-dependent (Yang et al., 2013). We aim to minimize the
reliance on a specific cosmological model. Furthermore, the uncertainty in the
luminosity distance of GRB data necessitates the use of a specific likelihood
function (Reichart et al., 2001), which is not pertinent to this study.

We use the `scikit-learn` 111https://scikit-learn.org/stable/index.html module
to demonstrate the general GP reconstruction generated using LML training
hyperparameters. The hyperparameters are optimized by LML, where the constant-
value and length-scale correspond to $\sigma_{f}^{2}$ and $l$, respectively.
Figure 1 shows that different kernel selections result in distinct curves
after reconstruction, but it is challenging to infer which performs better
from the graphs. To address this issue, we propose using the Approximate
Bayesian Computation (ABC) Rejection technique. The results of the GP
reconstruction of the observed data will also be needed when using ABC
Rejection, where the GP will be additionally programmed instead of using the
one in `scikit-learn`.

Table 1: 32 $H(z)$ Measurements from the Cosmic Chronometer Method

$z$ | $H(z)$ | $\sigma_{H}$ | References
---|---|---|---
0.09 | 69 | 12 | Jimenez et al. (2003)
0.17 | 83 | 8 |
0.27 | 77 | 14 |
0.4 | 95 | 17 |
0.9 | 117 | 23 |
1.3 | 168 | 17 | Simon et al. (2005)
1.43 | 177 | 18 |
1.53 | 140 | 14 |
1.75 | 202 | 40 |
0.48 | 97 | 62 | Stern et al. (2010)
0.88 | 90 | 40 |
0.1791 | 75 | 4 |
0.1993 | 75 | 5 |
0.3519 | 83 | 14 |
0.5929 | 104 | 13 |
0.6797 | 92 | 8 | Moresco et al. (2012)
0.7812 | 105 | 12 |
0.8754 | 125 | 17 |
1.037 | 154 | 20 |
0.07 | 69 | 19.6 |
0.12 | 68.6 | 26.2 | Zhang et al. (2014)
0.2 | 72.9 | 29.6 |
0.28 | 88.8 | 36.6 |
1.363 | 160 | 33.6 | Moresco (2015)
1.965 | 186.5 | 50.4 |
0.3802 | 83 | 13.5 |
0.4004 | 77 | 10.2 |
0.4247 | 87.1 | 11.2 | Moresco et al. (2016)
0.4497 | 92.8 | 12.9 |
0.4783 | 80.9 | 9 |
0.47 | 89 | 34 | Ratsimbazafy et al. (2017)
0.80 | 113.1 | 28.5 | Jiao et al. (2022)

Note. — The observational Hubble parameter $H(z)$ measurements of CC (in units
of km$\cdot$s${}^{-1}\cdot$Mpc-1) and their errors $\sigma_{H}$ at redshift
$z$. Here we have 32 data points in total.

Figure 1: Reconstructed images of GP obtained using scikit-learn. For all 3
graphs, the different coloured curves represent the curves of the different
kernel reconstructions, the red data points with error bars represent the
observed data and the shaded regions corresponding to the colour of the curves
indicate the range of 95% confidence intervals of the reconstruction. For each
of them, we have the following descriptions.
(a): Reconstruction for CC. The horizontal axis $z$ represents the redshift
and the vertical axis $H(z)$ represents the value of the Hubble parameter
obtained from the observation.
(b): Reconstruction for SNIa(binned). The horizontal axis $\ln(z)$ represents
the redshift value after logarithmization and the vertical axis $m(z)$
represents the distance modulus of the supernovae.
(c): Plot for reconstruction of GP using GRB dataset. The horizontal axis has
the same meaning as in Fig.(b). The vertical axis $\mu(z)$ represents the
distance modulus of the GRB event.

###  2.2 Approximate Bayesian Computation(ABC) Rejection

The foundation of the standard Bayesian inference procedure lies in the prior
distribution of the parameters and the likelihood function, which use the
Bayes rule to determine the posterior distribution. The Bayes rule is
expressed as follows:

$P(\theta|\boldsymbol{y})=\frac{P(\theta)P(\boldsymbol{y}|\theta)}{\int
P(\boldsymbol{y}|\theta)P(\theta)\mathrm{d}\theta},$ (13)

where $P(\theta|\boldsymbol{y})$ denotes the posterior distribution, which is
obtained through Bayesian analysis. Given the observed dataset
$\boldsymbol{y}$, it represents the probability distribution of the parameter
$\theta$ value. The researcher can determine the prior distribution of the
parameter $\theta$ denoted by the prior $P(\theta)$. The prior is the
probability function, or $P(\boldsymbol{y}|\theta)$. The normalizing constant
is $\int P(\boldsymbol{y}|\theta)P(\theta)\mathrm{d}\theta$, which is referred
to as “evidence”. For simplicity, we leave it at 1. Bayes rule is currently
expressed as:

$P(\theta|\boldsymbol{y})=P(\theta)P(\boldsymbol{y}|\theta).$ (14)

Bayesian analysis, based on Equation (14), when combined with the MCMC
approach, has been widely used in cosmological parameter estimation (CPE) and
has become a standard way of analysis in CPE(Zuntz et al., 2015). Bayesian
analysis can also be used for model selection(Trotta, 2008). Furthermore, in
our work, choosing kernels can be seen as a form of model selection, and thus,
our conception of model selection is also founded on Bayesian analysis. When
applying Bayesian analysis to model selection, Equation (14) must be modified
to take the following form:

$P(M_{i}|\boldsymbol{y})=P(M_{i})P(\boldsymbol{y}|M_{i}).$ (15)

The posterior distribution that we require represents the probability of the
collection of observations being produced by the model $M_{i}$. However,
directly applying Equation (15) is difficult since the likelihood function
$P(\boldsymbol{y}|M_{i})$ cannot be represented analytically when the various
kernels are viewed as distinct models.

It is important to note that although GP assumes that the reconstructed
functions in the model follow a Gaussian likelihood distribution, this does
not imply that $P(\boldsymbol{y}|M_{i})$ is a Gaussian likelihood. Here
$M_{i}$ represents different kernels. In fact, GP assumes that the vector of
observations is obtained from a Gaussian distribution in the function
space(Seikel et al., 2012), and that the mean of this Gaussian distribution is
our prediction $\overline{\boldsymbol{f^{*}}}$. Hence, if we use Bayes’
theorem for kernel selection, the resolved form of its likelihood function
remains unknown to us. This is also mentioned in Bernardo & Said (2021). Since
we cannot represent the likelihood function analytically, the traditional MCMC
approach is not applicable to our problem.

As an alternative, we choose Approximate Bayesian Computation (ABC) Rejection,
which has the benefit of not requiring the specification of a likelihood
function(Turner & Van Zandt, 2012). This quality makes the ABC Rejection
approach popular in cosmology (Jennings & Madigan (2017);Weyant et al.
(2013);Ishida et al. (2015)). In ABC Rejection, the likelihood function is
approximated by using frequencies to estimate probabilities, allowing us to
derive the posterior distribution. Here, we sample the model’s parameters
multiple times, referring to each sample as a particle. We then establish
suitable screening criteria and calculate the fraction of particles that pass
the screening relative to the total number of samples to derive the frequency
and, consequently, the probability. To apply the ABC Rejection, we treat the
kernel function as a model $M$, and sample the hyperparameters
$\sigma_{f}^{2}$ and $l$ as parameters in the model $M$. The following are the
steps to use the ABC Rejection method(Toni & Stumpf, 2009):

  * •

To ensure uniformity across all kernels studied, we assign the same prior
probability, $p(M_{i})$, to each kernel. From this collection of kernels, we
select a specific kernel, $k(z_{i},z_{j})$, as our model $M$. The
hyperparameters, $(\sigma_{f}^{2},l)$, which are treated as a ’particle’, are
then sampled from the prior distribution.

  * •

Using the sampled hyperparameters $(\sigma_{f}^{2},l)$ and the selected kernel
$k(z_{i},z_{j})$, we obtain the mean of the reconstructed function, denoted by
$\overline{\boldsymbol{f^{*}}}$. This reconstructed function is generated
using GP. It is important to note that the mean of the reconstructed function
$\overline{\boldsymbol{f^{*}}}$ is equivalent to a simulated dataset.

  * •

We then calculate the distance, denoted by $d$, between the simulated dataset
$\overline{\boldsymbol{f^{*}}}$ and the observed dataset $\boldsymbol{y}$ by
substituting both into the distance function $d(x,y)$.

  * •

To determine whether the sampled hyperparameters $(\sigma_{f}^{2},l)$ are
acceptable, we set a threshold, denoted by $\varepsilon$, and compare the
calculated distance $d$ with $\varepsilon$. If $d\leq\varepsilon$, we accept
$(\sigma_{f}^{2},l)$ and add it to the count. Otherwise, we discard
$(\sigma_{f}^{2},l)$ and do not count it.

  * •

By repeating the above steps, we can obtain the posterior distribution
$p(\boldsymbol{y}|M_{i})$ that we want to calculate. Specifically, we sample
$(\sigma_{f}^{2},l)$ multiple times from a particular distribution and count
the number of particles that satisfy $d\leq\varepsilon$. Finally, we solve for
the ratio of that number to the total number of samples. When $\varepsilon$ is
sufficiently small to indicate obvious differences, this ratio corresponds to
the posterior distribution $p(\boldsymbol{y}|M_{i})$.

As can be observed from the above stages, when computing the posterior
distribution using ABC Rejection for the results, the selection of the
distance function $d(x,y)$, the distribution of the hyperparameters sampled,
and the sampling technique chosen are all critical. In our research, we
utilized hyperparameters from a uniform distribution, which is available in
the `abcpmc` package(Akeret et al., 2015).

In practice, many previous studies combined importance sampling algorithms to
assign weights to the particles, gradually reducing the threshold to conserve
computational resources for the ABC Rejection algorithm. The SMC (Sequential
Monte Carlo) algorithm(Toni & Stumpf, 2010) and the PMC (Population Monte
Carlo) algorithm(Cappé et al., 2004) are two typical algorithms that are
computed iteratively, and the threshold is reduced after each generation
according to the weight of the particles that have not yet been filtered out.
Therefore, adaptive algorithms for weight calculation can also be
applied(Beaumont et al. (2009);Bonassi & West (2015)). Although the ABC-PMC
algorithm is usually employed for parameter inference without likelihood
functions(Kilbinger et al., 2012) and less frequently for model selection, the
ABC-SMC algorithm is very effective when dealing with large-scale data and has
also been used for model selection(Abdessalem et al. (2017);Bernardo & Said
(2021)).

To assign weights to the particles, several previous studies integrated
importance sampling approaches, gradually lowering the threshold to conserve
computational resources for the ABC Rejection procedure. The SMC (Sequential
Monte Carlo) algorithm(Toni & Stumpf, 2010) and the PMC (Population Monte
Carlo) algorithm(Cappé et al., 2004) are two examples of these algorithms.
They are calculated iteratively, and the threshold is lowered in accordance
with the weight of the particles that haven’t yet been filtered out after each
generation. As a result, adaptive algorithms can also be used to calculate
weights (see Beaumont et al. (2009);Bonassi & West (2015)). While the ABC-PMC
algorithm is typically employed for parameter inference without likelihood
functions(Kilbinger et al., 2012) and less frequently for model selection, the
ABC-SMC method has been used for model selection and is particularly efficient
when working with large datasets. However, since the dataset we use is not
particularly large, and the ABC-SMC technique only returns ”which is the best
kernel” and cannot compare the magnitude of each kernel relative to the data,
we do not intend to employ it in our study. In fact, even without importance
sampling in our investigation, we can still obtain consistent results. We only
need to consider how different distance functions affect the posterior
distribution. Section 2.3 describes the distance functions used in the study.

###  2.3 Distance Function

The selection of an appropriate distance function is critical for Approximate
Bayesian Computation (ABC) analysis as different functions may indicate
varying levels of statistical significance between simulated and real
datasets. Previous research in this field has typically used the likelihood
function (LML) as the distance function (Abdessalem et al. (2017); Bernardo &
Said (2021)). LML is frequently used to evaluate the impact of hyperparameter
values on the fit, making it a reasonable choice as one of the distance
functions. However, selecting hyperparameter values is not limited to LML
optimization and other methods can also be used (Gómez-Valent & Amendola,
2018). In our work, we selected several distance functions, including LML,
$\chi^{2}$ estimation, and bias, to compare their performance in ABC analysis:

  * •

$\chi^{2}$(Bernardo & Said, 2021)

$\displaystyle\chi^{2}=d(\overline{\boldsymbol{f^{*}}},\boldsymbol{y})=(\overline{\boldsymbol{f^{*}}}-\boldsymbol{y})^{T}\Sigma_{y}^{-1}(\overline{\boldsymbol{f^{*}}}-\boldsymbol{y}),$
(16)

  * •

Bias

$\displaystyle
d(\overline{f^{*}},y)=\frac{\sum_{i=1}^{n}|y_{i}-\overline{f^{*}}_{i}|}{n}.$
(17)

Optimizing the residuals of the data is equivalent to minimizing Equation
(16). This equation considers the need to minimize the sum of squares of the
residuals while taking into account the weighting of the inverse error. Based
on the residuals, Equation (16) provides a criterion for assessing the model’s
merit, where a smaller value of $\chi^{2}$ indicates better agreement between
the simulated and observed data.

Equation (17) provides the average of Euclidean distances between the
simulated and observed datasets and represents the expectation of the
difference between the projected and actual values of the model, also known as
bias. The value of bias indicates how well the model fits the data, and the
lower its value, the closer the mean value of the simulated data is to the
observed data.

By combining these three distance functions, we provide three different
methods for filtering particles, and the ABC Rejection results obtained from
these methods offer a more comprehensive answer to the question of which
kernel performs best in ABC analysis.

###  2.4 Bayes Factor

The Bayes factor, denoted by $K$, is a measure of the support for one
statistical model over another, obtained by taking the ratio of their
likelihoods(Morey et al., 2016). In our work, we use the Bayes factor to
assess the dependence of different datasets on the kernel. Unlike traditional
hypothesis testing, which only allows for acceptance or rejection of a
hypothesis, the Bayes factor evaluates the evidence in support of a
hypothesis. Hence, the Bayes factor not only identifies the best model among
competing kernels but also quantifies how much better it is than the others.
For two alternative models $M_{1}$ and $M_{2}$ with observed data
$\boldsymbol{y}$, the plausibility of the models is evaluated by the Bayes
factor $K$, which is given by Equation (18):

$K_{12}=\frac{P(\boldsymbol{y}|M_{1})}{P(\boldsymbol{y}|M_{2})}=\frac{P(M_{1}|\boldsymbol{y})}{P(M_{2}|\boldsymbol{y})}\frac{P(M_{2})}{P(M_{1})}.$
(18)

In fact, the Bayes factor is the ratio of the evidence from the two models.
But in order to be able to use the results of ABC Rejection, we have
transformed the original form of the Bayes factor into Equation (18) to use
the results of ABC Rejection (see Section 2.5 for a direct calculation of the
evidence). When the prior probabilities of the two models are equal
($P(M_{1})=P(M_{2})$), the Bayes factor $K$ represents the ratio of the
posterior probability of $M_{1}$ and $M_{2}$. We compute the same prior
probability for the two kernels when we calculate the Bayes factor between
them. The procedure only takes into account the ratio of the posterior
distributions of the two kernels.

In Equation (18), a value of $K>1$ implies that $M_{1}$ is more supported by
the data than $M_{2}$. The scale of $K$ has a quantitative interpretation
based on probability theory(Jeffreys, 1998), as shown in Table 2. We use this
table to determine the strength of evidence for each of the two kernels.

Table 2: Correspondence between the Value of $K$ and the Strength of the
Evidence

$K_{12}$ | $\log_{10}K$ | Strength of Evidence
---|---|---
$<10^{0}$ |  $<$0 | Negative (supports $M_{2}$)
$10^{0}$ to $10^{1/2}(\approx 3.162)$ | 0 to 1/2 | Barely worth mentioning
$10^{1/2}(\approx 3.162)$ to $10^{1}$ | 1/2 to 1 | Substantial
$10^{1}$ to $10^{3/2}(\approx 31.623)$ | 1 to 3/2 | Strong
$10^{3/2}(\approx 31.623)$ to $10^{2}$ | 3/2 to 2 | Very strong
$>10^{2}$ |  $>$2 | Decisive

Note. — $K$ denotes the value of the Bayes factor and $\log_{10}K$ denotes the
result of taking the usual logarithm of $K$. The third column indicates the
strength of the evidence in the corresponding range of $K$ values. Subsequent
analyses of the Bayes factor results are based on this table. This table gives
the criteria for quantitatively evaluating the relative reliability between
two models.

###  2.5 Nested Sampling

Nested sampling is a method originally proposed by John Skilling to
approximate the evidence in Bayes’ theorem(Skilling, 2006). Bayes’ theorem
takes the form of:

$\displaystyle P(\theta|D,M)$
$\displaystyle=\frac{P(D|\theta,M)P(\theta|M)}{P(D|M)},$ $\displaystyle
P(D|M)$ $\displaystyle=\mathcal{Z}=\int
P(D|\theta,M)P(\theta|M)\mathrm{d}\theta.$ (19)

In Bayes’ theorem, $P(D|M)$ is referred to as the “evidence” and is an
unresolvable representation of the likelihood function, as mentioned in ABC
Rejection. $D$ represents our observed data, while $\theta$ represents the
parameter of the model $M$ (which may include additional parameters). To make
model comparisons, the easiest approach is to compare the size of the evidence
corresponding to different models with the same prior probability. Nested
sampling is the method by which we can obtain the evidence directly.

To perform model comparisons using nested sampling, we first rewrite Bayes’
theorem as follows(Handley et al., 2015a):

$\displaystyle P(M_{i}|D)$ $\displaystyle=\frac{P(D|M_{i})P(M_{i})}{P(D)},$
$\displaystyle=\frac{\mathcal{Z}_{i}\pi_{i}}{\Sigma_{j}\mathcal{Z}_{j}\pi_{j}}.$
(20)

In Equation (2.5), $\pi$ refers to the priors of the model $M_{i}$. In
practice, it is often assumed that $\pi$ is equal for each model to be
compared, and so Equation (2.5) is not necessary for model comparison. The
only integral that needs to be calculated is that expressed in Equation (2.5),
and the model with the greater value of evidence is deemed to be superior in
the comparison.

The principle of nested sampling can be briefly described as follows(Handley
et al., 2015a). We begin by sampling $n_{\text{live}}$ points from the prior
distribution. We then define a prior volume $X(\mathcal{L})$:

$X(\mathcal{L})=\int_{\mathcal{L}(\theta)>\mathcal{L}}\pi(\theta)\mathrm{d}\theta,$
(21)

where $\mathcal{L}$ is a likelihood we set. At each iteration $i$, we remove
the points with likelihood lower than $\mathcal{L}$ and replace them with new
live points. According to Equation (2.5), the integral now takes the form of:

$\mathcal{Z}=\int_{0}^{1}\mathcal{L}(X)\mathrm{d}X.$ (22)

The integral can be written in the form of a summation to the point to get the
approximation of the evidence. More details are shown in Skilling (2006),
Handley et al. (2015a) and Handley et al. (2015b).

To implement nested sampling, we can use the Python package `Dynesty`(Speagle,
2020), which is easy to install and use. With `Dynesty`, we input a prior and
a log-marginal-likelihood, and its internal nested sampler is called to
complete the sampling. For GP, the prior used is not the same as the prior of
hyperparameters; it refers to the reconstructed data $\boldsymbol{y}$.
However, the log-marginal-likelihood is the same. In `Dynesty`, the log-
marginal-likelihood takes the form of $P(D|\theta,M)$, and for GP, the LML is
$P(\boldsymbol{y}|\sigma_{f},l,\boldsymbol{X})$. Here, $\boldsymbol{X}$ is
equal to $\theta$, while $\sigma_{f}$ and $l$ describe the model $M$. Thus, we
can directly use the LML in GP as the likelihood in `Dynesty`.

Another important option in `Dynesty` is the choice of a specific bounding. We
choose “Multiple ellipsoid bounds”, as it is the most frequently used
method(Feroz et al., 2009).

In `Dynesty`, the program output for evidence is usually $\ln\mathcal{Z}$. As
the essence of the Bayes factor is the division of the evidence from two
models, $\ln\mathcal{Z}$ becomes a subtraction between different models’
$\ln\mathcal{Z}$:

$\ln K_{12}=\ln P(\boldsymbol{y}|M_{1})-\ln
P(\boldsymbol{y}|M_{2})=\ln\mathcal{Z}_{1}-\ln\mathcal{Z}_{2}.$ (23)

As in Table 2, such a comparison table will be given for $\ln K_{12}$. The
relationship between this difference and the strength of model support is
shown in the Tabel 3(Sarro et al., 2012).

Table 3: Correspondence between the Value of $\ln K_{12}$ and the Strength of the Evidence $|\ln K_{12}|$ | odds | Strength of evidence
---|---|---
$<1.0$ | $\lesssim 3:1$ | Inconclusive
$1.0$ | $\sim 3:1$ | Weak evidence
$2.5$ | $\sim 12:1$ | Moderate evidence
$5.0$ | $\sim 150:1$ | Strong evidence

##  3 Results

###  3.1 ABC Rejection with Different Distance Functions

In this investigation, we employ three different types of data, apply the ABC
rejection method to each type of data, and use three different distance
functions in the computations, resulting in nine graphs presented in Figure 2.

Figure 2: For each subplot, we utilize a bar plot format to display the
posterior distribution of the five kernels. The horizontal axis for each
subplot denotes the type of kernel, while the vertical axis indicates the
likelihood that the observed data is created by them. Normalization is applied
to the probability.

It is important to note that the size of the threshold is not chosen
arbitrarily. Setting the threshold too high would obscure the differences
between specific kernel functions, while setting the threshold too low would
result in not only a small number of particles in each kernel function, but
also particles that are very close to one another. For instance, when the
threshold is set too low, the RBF kernel produces only two particles that are
closer to the observed data than the threshold distance, while the M52 kernel
produces four particles. Although the posterior distribution of the M52 kernel
is twice as large as that of the RBF kernel in terms of normalized
probability, the difference between the two is not significant. Therefore, an
incorrect judgment would be made if the threshold is set too low.

To address this issue, we continuously adjust the threshold until we reach the
final result. When the posterior distributions of the individual kernels
undergo significant changes when the threshold is set to $\varepsilon$, but do
not differ significantly when the threshold is greater than $\varepsilon$, we
consider $\varepsilon$ to be the appropriate threshold. When the previously
observed differences are preserved when the threshold is set to a value less
than $\varepsilon$, we consider $\varepsilon$ to be the correct threshold. It
is worth mentioning that there are no circumstances where the differences in
the posterior distributions of the kernels change when the threshold is
decreased further.

After resolving the threshold selection problem, our next step is to reduce as
much randomness as possible from the results. We ensure that the sampling
opportunities are equal for all 5 kernels, meaning that their prior
probabilities are all 20%. This is reflected in the sampling of $2\times
10^{5}$ samples for each kernel in each scenario. At the end of every $10^{4}$
samples, we record each kernel’s posterior distribution, resulting in 100
posterior probabilities for each kernel in each example. The posterior
distribution for each kernel function in Figure 2 is obtained by averaging
these 100 posterior probabilities.

As shown in Figure 2, the M52 kernel has the largest posterior distribution in
all situations, indicating that the GP reconstructed by the M52 kernel is the
closest to the real data and is therefore the best model. The extent to which
M52 outperforms the other kernels varies between subplots.

It is worth noting that the RBF kernel, which is commonly used by astronomers,
has a larger bias in all reconstruction instances than the other kernels. On
the other hand, the M52 kernel reconstruction fits the data much better than
the RBF kernel.

Based on the ABC Rejection results, we can conclude that M52 is the best
kernel for the purpose of this study. Similar results are presented in
Bernardo & Said (2021), but there are no results for M52 compared to other
models. The strength of the evidence is discussed in Section 3.2 in comparison
to other kernels.

###  3.2 Bayes Factor $K$ Between Different Kernels (in ABC Rejection)

We generate heat maps to visualize the magnitude of the Bayes factor $K$
between the kernels for each of the nine cases depicted in Figure 2. A
colorbar is used to represent the Bayes factor values, with significant values
of 1, 3.16, and 10 labeled.

Figure 3: The Bayes factor between every 2 kernels is visualized using this
heat map to represent the strength of evidence level of the observed data
between any two kernels. As an example, the value corresponding to the color
block in the second column of the first row on the colorbar indicates the
Bayes factor of CHY with respect to RBF. That is,
$P(M_{\text{CHY}}|\boldsymbol{y})/P(M_{\text{RBF}}|\boldsymbol{y})$(According
to Equation (18)). The dataset and distance function related to the posterior
distribution utilized for the Bayes factor computation for that subplot are
indicated in the subplot’s title.

Figure 3 shows that the M52 kernel significantly outperforms the RBF kernel,
while the difference in performance between M52 and the other kernels is not
significant, except for the comparison with RBF.

For the CC dataset, the strength of evidence for M52 is only at the “Barely
worth mentioning” level compared to RBF and other kernels. However, the Bayes
factor of M52 for RBF is the highest among all cases.

When using the SNIa dataset, the evidence strength of M52 is at the “Strong”
level compared to RBF, occasionally reaching “Strong” when compared to other
kernels but mostly staying at “Substantial” and “Barely worth mentioning”. The
M72 kernel is the second-best performing kernel after M52.

For the GRB dataset, the evidence strength for M52 relative to RBF is
“Substantial” with LML, Bias, and $\chi^{2}$ distance functions. Other kernels
have evidence ranging from “barely worth mentioning” to “strong”, with M72
being the second-best performing kernel after M52.

Although the M52 kernel outperforms all other kernels, the strength of the
evidence is not always decisive. Therefore, in most cases, replacing M52 with
a different kernel would not affect the GP reconstruction significantly.
However, when choosing between M52 and RBF, we prefer to use M52 as the
commonly used kernel in GP.

###  3.3 Nested Sampling: Evidence and Bayes factors

The nested sampler of `Dynesty` is utilized in this study to produce Figure 4.
For each type of data, four figures are generated by the sampler. The bottom
graph of each column in the figure is particularly useful for model selection,
as it displays how the evidence for the model changes over the sampling
process until it stabilizes. The other three plots in each column have
different meanings, which are more thoroughly explained in Speagle (2020) and
are not discussed in detail in our work.

Figure 4: Plots for sampler’s running of 3 datasets. The horizontal axis of
each subplot represents the natural logarithm of the a prior volume, $-\ln X$.
The plot with the vertical axis of $\ln\mathcal{Z}$ is the plot we need for
model comparison, and the value that tends to be stable is the model’s
corresponding evidence $\ln\mathcal{Z}$.

Remarkably, our results differ from those of ABC Rejection. In Figure 4, it is
evident that M52 no longer dominates and even exhibits the lowest evidence for
the CC data. The SNIa data, on the other hand, does not present an obvious
result from the graph. To compare the results, we extract the exact runs’
evidence results and display them in Table 4. We then visualize these results
using the bar chart in Figure 5.

Table 4: The evidence $\ln\mathcal{Z}$ obtained for all types of data and kernels Data | Kernels | $\ln\mathcal{Z}$
---|---|---
| RBF | -128.839$\pm$0.053
| CHY | -129.015$\pm$0.057
CC | M52 | -129.494$\pm$0.060
| M72 | -128.875$\pm$0.052
| M92 | -128.764$\pm$0.051
| RBF | 85.729$\pm$0.705
| CHY | 85.171$\pm$0.863
SNIa | M52 | 85.882$\pm$0.751
| M72 | 84.745$\pm$0.843
| M92 | 85.746$\pm$0.762
| RBF | -164.998$\pm$0.140
| CHY | -164.631$\pm$0.138
GRB | M52 | -165.449$\pm$0.149
| M72 | -165.166$\pm$0.145
| M92 | -164.703$\pm$0.140
Figure 5: Visualization plots of evidence in all cases. The solid black line
indicates the error bar. The errors are consistent with those given in Table
4. Note that for the CC and GRB data, the log evidence for them is negative.

From Figure 5, it is apparent that, except for the SNIa data where M52 is the
best kernel, the other two cases show the opposite result, i.e. M52 is the
“worst” kernel. Even for the SNIa data, M52 does not display a decisive
advantage over the other kernels.

It is also important to understand how the nested sampling results reflect the
Bayes factor’s strength of evidence between different kernels, as discussed in
Equation (23). To display the Bayes factors between the kernels, similar to
the approach used in ABC Rejection (Figure 3), we also generated a heatmap in
Figure 6.

To make the Bayes factor size associated with each color block identifiable,
we label them in Figure 6. However, since the difference between each kernel
is minor, the Bayes factors do not display any significant differences. The
largest Bayes factor only reaches the level of “weak evidence” (M92 vs. M72
for SNIa data). Nearly all Bayes factors remain at the “inconclusive” level,
as shown in Table 3. Therefore, we can assume that there is no significant
difference among commonly used kernels in GP. Researchers can select kernel
functions for their practical studies without concern for significant
deviations in the results.

Figure 6: Bayes factor heat map obtained by direct computation of the
evidence. The meaning of each color block is the same as in Figure 3. As an
example, for the color block in the second row and first column, the Bayes
factor it represents is obtained from $\ln\mathcal{Z}_{\text{RBF}}-$
$\ln\mathcal{Z}_{\text{CHY}}$.

##  4 conclusion

In this study, we apply ABC Rejection to investigate the topic of kernel
selection for Gaussian processes (GP) using three alternative distance
functions. These distance functions correspond to different ways of particle
filtering using the ABC approach. Our results demonstrate that even without
the additional sampling method of lowered thresholding, reliable kernel
selection outcomes may be obtained. We find that the M52 kernel performs as
the optimal kernel in all nine scenarios studied. Additionally, the widely
employed RBF kernel function by astronomers does not perform as well as the
other kernels, especially when compared to M52. A recent paper (Bernardo &
Said, 2021) reached a similar conclusion, but their analysis was limited to CC
and SNIa data only and did not include a comparison of M52 to other kernels.
Our research is more comprehensive and thorough, yielding more definitive
results.

Moreover, the Bayes factor calculation results indicate that M52 outperforms
the other options. The results reveal that, according to Table 2, M52 can only
attain a “strong” level of evidence for the observed data in a few
circumstances when compared to the RBF kernel. However, M52 does not exhibit
very strong support for the observed data when compared to the other options.
As a result, while M52 is the best kernel function in our research range,
there is little statistical difference in the outputs produced by using
kernels other than the RBF kernel.

On the other hand, we use a nested sampling algorithm and the `Dynesty`
package to directly evaluate the evidence for each kernel function to compare
their performance under different data. This approach leads to different
results from the ABC approach, where M52 no longer maintains its dominance in
the evidence comparison, but instead becomes the kernel function with the
lowest evidence in the CC and GRB data. However, this does not significantly
affect our conclusions. By taking the log Bayes factors between them, we find
that these differences are not remarkable, and no kernel function is
appreciably supported by the data. Thus, based on the results of nested
sampling, it does not matter which of these commonly used kernel functions a
researcher chooses when using GP (although this conclusion is specific to the
type of dataset we investigated).

There are several other interesting points to consider regarding nested
sampling. First, the computational speed of the program is a point of interest
for improvement. The slow speed of the program using the `Dynesty` package
makes it challenging to run the algorithm multiple times to rule out
randomness in the results. In our work, the ABC Rejection algorithm runs a
million times for each kernel, and the final result is obtained by averaging
over bins (100 times for a bin). However, repeated runs of this size cannot be
performed for nested sampling. Second, we find that although repeated nested
sampling leads to small fluctuations in the evidence for each kernel function,
for both CC and GRB data, M52 is always the kernel function with the lowest
evidence, and the results vary from time to time as to which one has the
highest evidence. In contrast, algorithms such as ABC Rejection and similar
ABC-SMC (Bernardo & Said, 2021) give the exact opposite results. Although this
does not affect the final selection, this result, which is unlikely to be the
result of statistical error, deserves further investigation.

As a future direction of this study, there are potential improvements to be
made in the approach of finding posterior distributions. In addition to ABC
Rejection, the posterior distribution can also be generated using VBMC(Acerbi,
2018). This is achieved by progressively modifying a given distribution to
resemble the true distribution based on the observed data. To further limit
the selection of kernel, the introduction of more simulated data to populate
the observed data is possible(Bernstein et al., 2012). Another approach that
could be considered is the use of Sparse Markovian Gaussian
processes(Wilkinson et al., 2021). This method utilizes several local Gaussian
terms to approximate a non-Gaussian likelihood function, making it
particularly suitable for large sample data. As the likelihood function of
cosmological data is also uncertain, the issue of kernel selection for Sparse
Markovian Gaussian processes is of great significance.

We express our gratitude to the anonymous reviewer for the constructive
comments, and to Bo Yu, Xu-Yao Gong, Jing Niu, Yuan-Hang Mao, Yuan-Bo Xie, and
Kang Jiao for their insightful discussions. We would also like to acknowledge
the contributions of the editors. The python packages `NumPy`(Harris et al.,
2020) and `Matplotlib`(Hunter, 2007) played a critical role in our research,
and we extend our thanks to their developers. This work is supported by grants
from the National Science Foundation of China (Nos.11929301 and 61802428)

## References

  * Abdessalem et al. (2017) Abdessalem, A. B., Dervilis, N., Wagg, D. J., & Worden, K. 2017, Frontiers in Built Environment, 3, 52
  * Acerbi (2018) Acerbi, L. 2018, Advances in Neural Information Processing Systems, 31
  * Akeret et al. (2015) Akeret, J., Refregier, A., Amara, A., Seehars, S., & Hasner, C. 2015, Journal of Cosmology and Astroparticle Physics, 2015, 043
  * Beaumont et al. (2009) Beaumont, M. A., Cornuet, J.-M., Marin, J.-M., & Robert, C. P. 2009, Biometrika, 96, 983
  * Bernardo & Said (2021) Bernardo, R. C., & Said, J. L. 2021, Journal of Cosmology and Astroparticle Physics, 2021, 027
  * Bernstein et al. (2012) Bernstein, J., Kessler, R., Kuhlmann, S., et al. 2012, The Astrophysical Journal, 753, 152
  * Bonassi & West (2015) Bonassi, F. V., & West, M. 2015, Bayesian Analysis, 10, 171
  * Cappé et al. (2004) Cappé, O., Guillin, A., Marin, J.-M., & Robert, C. P. 2004, Journal of Computational and Graphical Statistics, 13, 907
  * Demianski et al. (2017) Demianski, M., Piedipalumbo, E., Sawant, D., & Amati, L. 2017, Astronomy & Astrophysics, 598, A112
  * Dhawan et al. (2021) Dhawan, S., Alsing, J., & Vagnozzi, S. 2021, Monthly Notices of the Royal Astronomical Society: Letters, 506, L1
  * Feroz et al. (2009) Feroz, F., Hobson, M., & Bridges, M. 2009, Monthly Notices of the Royal Astronomical Society, 398, 1601
  * Gómez-Valent & Amendola (2018) Gómez-Valent, A., & Amendola, L. 2018, Journal of Cosmology and Astroparticle Physics, 2018, 051
  * Handley et al. (2015a) Handley, W., Hobson, M., & Lasenby, A. 2015a, Monthly Notices of the Royal Astronomical Society, 453, 4384
  * Handley et al. (2015b) —. 2015b, Monthly Notices of the Royal Astronomical Society: Letters, 450, L61
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Hwang et al. (2023) Hwang, S.-g., L’Huillier, B., Keeley, R. E., Jee, M. J., & Shafieloo, A. 2023, Journal of Cosmology and Astroparticle Physics, 2023, 014
  * Ishida et al. (2015) Ishida, E. E., Vitenti, S. D., Penna-Lima, M., et al. 2015, Astronomy and Computing, 13, 1
  * Jeffreys (1998) Jeffreys, H. 1998, The theory of probability (OUP Oxford)
  * Jennings & Madigan (2017) Jennings, E., & Madigan, M. 2017, Astronomy and computing, 19, 16
  * Jiao et al. (2022) Jiao, K., Borghi, N., Moresco, M., & Zhang, T.-J. 2022, arXiv preprint arXiv:2205.05701
  * Jimenez et al. (2003) Jimenez, R., Verde, L., Treu, T., & Stern, D. 2003, The Astrophysical Journal, 593, 622, doi: 10.1086/376595
  * Kilbinger et al. (2012) Kilbinger, M., Benabed, K., Cappé, O., et al. 2012, Astrophysics Source Code Library, ascl
  * Mehrabi & Basilakos (2020) Mehrabi, A., & Basilakos, S. 2020, The European Physical Journal C, 80, 1
  * Moresco (2015) Moresco, M. 2015, Monthly Notices of the Royal Astronomical Society: Letters, 450, L16, doi: 10.1093/mnrasl/slv037
  * Moresco et al. (2012) Moresco, M., Cimatti, A., Jimenez, R., et al. 2012, Journal of Cosmology and Astroparticle Physics, 2012, 006, doi: 10.1088/1475-7516/2012/08/006
  * Moresco et al. (2016) Moresco, M., Pozzetti, L., Cimatti, A., et al. 2016, Journal of Cosmology and Astroparticle Physics, 2016, 014, doi: 10.1088/1475-7516/2016/05/014
  * Morey et al. (2016) Morey, R. D., Romeijn, J.-W., & Rouder, J. N. 2016, Journal of Mathematical Psychology, 72, 6, doi: https://doi.org/10.1016/j.jmp.2015.11.001
  * Mukherjee & Banerjee (2021) Mukherjee, P., & Banerjee, N. 2021, The European Physical Journal C, 81, 1
  * Ratsimbazafy et al. (2017) Ratsimbazafy, A. L., Loubser, S. I., Crawford, S. M., et al. 2017, Monthly Notices of the Royal Astronomical Society, 467, 3239, doi: 10.1093/mnras/stx301
  * Reichart et al. (2001) Reichart, D. E., Lamb, D. Q., Fenimore, E. E., et al. 2001, The Astrophysical Journal, 552, 57
  * Sarro et al. (2012) Sarro, L. M., Eyer, L., O’Mullane, W., & De Ridder, J. 2012, Astrostatistics and Data Mining, Vol. 2 (Springer Science & Business Media)
  * Scolnic et al. (2018) Scolnic, D. M., Jones, D., Rest, A., et al. 2018, The Astrophysical Journal, 859, 101
  * Seikel & Clarkson (2013) Seikel, M., & Clarkson, C. 2013, arXiv preprint arXiv:1311.6678
  * Seikel et al. (2012) Seikel, M., Clarkson, C., & Smith, M. 2012, Journal of Cosmology and Astroparticle Physics, 2012, 036
  * Shafieloo et al. (2012) Shafieloo, A., Kim, A. G., & Linder, E. V. 2012, Phys. Rev. D, 85, 123530, doi: 10.1103/PhysRevD.85.123530
  * Simon et al. (2005) Simon, J., Verde, L., & Jimenez, R. 2005, Phys. Rev. D, 71, 123001, doi: 10.1103/PhysRevD.71.123001
  * Skilling (2006) Skilling, J. 2006, Bayesian Analysis, 1, 833 , doi: 10.1214/06-BA127
  * Speagle (2020) Speagle, J. S. 2020, Monthly Notices of the Royal Astronomical Society, 493, 3132
  * Stern et al. (2010) Stern, D., Jimenez, R., Verde, L., Kamionkowski, M., & Stanford, S. A. 2010, Journal of Cosmology and Astroparticle Physics, 2010, 008, doi: 10.1088/1475-7516/2010/02/008
  * Sun et al. (2021) Sun, W., Jiao, K., & Zhang, T.-J. 2021, The Astrophysical Journal, 915, 123
  * Toni & Stumpf (2009) Toni, T., & Stumpf, M. P. 2009, arXiv preprint arXiv:0910.4472
  * Toni & Stumpf (2010) —. 2010, Bioinformatics, 26, 104
  * Trotta (2008) Trotta, R. 2008, Contemporary Physics, 49, 71
  * Turner & Van Zandt (2012) Turner, B. M., & Van Zandt, T. 2012, Journal of Mathematical Psychology, 56, 69
  * Weyant et al. (2013) Weyant, A., Schafer, C., & Wood-Vasey, W. M. 2013, The Astrophysical Journal, 764, 116
  * Wilkinson et al. (2021) Wilkinson, W., Solin, A., & Adam, V. 2021, in International Conference on Artificial Intelligence and Statistics, PMLR, 1747–1755
  * Williams & Rasmussen (2006) Williams, C. K., & Rasmussen, C. E. 2006, Gaussian processes for machine learning, Vol. 2 (MIT press Cambridge, MA)
  * Yang et al. (2013) Yang, X., Yu, H.-R., Zhang, Z.-S., & Zhang, T.-J. 2013, The Astrophysical Journal Letters, 777, L24
  * Zhang et al. (2014) Zhang, C., Zhang, H., Yuan, S., et al. 2014, Research in Astronomy and Astrophysics, 14, 1221, doi: 10.1088/1674-4527/14/10/002
  * Zuntz et al. (2015) Zuntz, J., Paterno, M., Jennings, E., et al. 2015, Astronomy and Computing, 12, 45
  * Ó Colgáin & Sheikh-Jabbari (2021) Ó Colgáin, E., & Sheikh-Jabbari, M. M. 2021, The European Physical Journal C, 81, 892, doi: 10.1140/epjc/s10052-021-09708-2


# PA-GM: Position-Aware Learning of Embedding Networks for Deep Graph Matching

###### Abstract

Graph matching can be formalized as a combinatorial optimization problem,
where there are corresponding relationships between pairs of nodes that can be
represented as edges. This problem becomes challenging when there are
potential ambiguities present due to nodes and edges with high similarity, and
there is a need to find accurate results for similar content matching. In this
paper, we introduce a novel end-to-end neural network that can map the linear
assignment problem into a high-dimensional space augmented with node-level
relative position information, which is crucial for improving the method’s
performance for similar content matching. Our model constructs the anchor set
for the relative position of nodes and then aggregates the feature information
of the target node and each anchor node based on a measure of relative
position. It then learns the node feature representation by integrating the
topological structure and the relative position information, thus realizing
the linear assignment between the two graphs. To verify the effectiveness and
generalizability of our method, we conduct graph matching experiments,
including cross-category matching, on different real-world datasets.
Comparisons with different baselines demonstrate the superiority of our
method. Our source code is available under https://github.com/anonymous.

Index Terms—  Graph Matching, Graph Embedding, Deep Neural Network

##  1 Introduction

Graph matching aims to establish the relationship between two or more graphs
based on information derived from their nodes and edges [1]. Due to its
ability to better express and encode data relationships, graph matching has
gained considerable traction in computer vision and related fields.
Applications of graph matching techniques include image registration in
medical image analysis [2], link analysis in social networks [3] and image
extrapolation in computer vision[4]. Current methods for solving the graph
matching problem can be divided into two main classes: a) learning-free or b)
learning-based [5, 6]. In the first case, the critical element is the
mathematical method used for obtaining an approximate solution to an
intrinsically NP-hard problem. On the other hand, learning-based methods aim
to improve the solver with a number of methods including deep neural networks.

However, these methods can only compute a similarity score for the whole
graph, or rely on inefficient global matching procedures[7]. For example, [8]
only considers the embedding of local information for nodes in the graph,
leading to a tendency to inconsistently match similar nodes from different
regions of the graph, and thus resulted to have ambiguities. As is illustrated
in Fig. 1, it shows that the node embedding representation, which relies only
on local structural information and semantic node information, lacks
sufficient discrimination to effectively resolve these ambiguities. As a
result, it is difficult to distinguish the left and right ears of the dog in
the example. From this, relative positional information is the key to the
matching of diagrams, especially in cases where the graphs have similar
semantic and structural content.

Fig. 1: A failed example of mismatching the position of the animal dog’s ears
in the two images.

In this paper, we introduce the idea of position-awareness to solve the above-
mentioned problem in graph matching. We first construct the graph as input to
the model with node features extracted from an image, and then extract a
collection of anchor points as reference coordinates for each node in the
graph. We further propose a corresponding position-aware node embedding
algorithm to capture the relative positional information of the nodes in the
graph. With the node-wise graph embedding to hand, we identify a node
permutation for node-to-node correspondence.

Fig. 2: Overview of the end-to-end position embedding networks for deep graph
matching. The blue source graph $\mathcal{G}_{s}$ and the green target graph
$\mathcal{G}_{t}$ are extracted with the node-wise graph feature
representation in high-level through the two frameworks of position-aware node
embedding.

The main contributions of this paper are as follows:

  * •

We propose a novel Position-Aware Learning of Embedding Network for Deep Graph
Matching (PA-GM). To the best of our knowledge, there are no methods that
consider the learning of an embedding augmented with relative positional
information in graph matching tasks. This restriction hampers their
applications in terms of matching accuracy.

  * •

We extract the relative positional information for the nodes in constructing
the node embedding instead of the traditional neighborhood aggregation
approach. We fully consider the positional information relevant to the
matching of keypoints, which is proved to be effective in our experiments.

  * •

The proposed framework is both scalable and flexible, and benefits from the
use of a relative position coefficient together with an alignment loss.
Experiments show that the model achieves state-of-the-art results on real-
world datasets.

##  2 Methods

In this paper, we intend to resolve the graph matching problem based on the
supervised matching of graphs. Specifically, we aim to learn an end-to-end
model which can extract graph information and their matches through given
pair-wise ground-truth correspondences for a set of graphs, and which can be
further generalized to unseen graph pairs. The overall pipeline for graph
matching using the position-aware embedding network is presented in Fig. 2.

###  2.1 Problem Definition

In this paper, we denote a graph by the triple $\mathcal{G}=(V,A,X)$ which
consists of a finite set of nodes $V$, an adjacency matrix $A$, and a set of
node attributes $X$ extracted from images using CNN-based models[9, 10]. We
construct a vector $\mathbf{v}\in\\{0,1\\}^{nm\times 1}$ to indicate the match
of vertices in the source graph $\mathcal{G}_{s}=(V_{s},A_{s},X_{s})$ and the
target graph $\mathcal{G}_{t}=(V_{t},A_{t},X_{t})$. The vector has elements
$\mathbf{v}_{i,j}=1$ if vertex $i\in V_{s}$ is matched to vertex $j\in V_{t}$,
and $\mathbf{v}_{i,j}=0$ if otherwise. It is worth noting that all the vertex
matches are subject to one-to-one mapping constraints $\sum_{j\in
V_{t}}\mathbf{v}_{i,j}=1\;\forall i\in V_{s}$ and $\sum_{i\in
V_{s}}\mathbf{v}_{i,j}\leq 1\;\forall j\in V_{t}$. Furthermore, we construct a
square symmetric positive matrix $M\in\mathbb{R}^{nm\times nm}$ as the
affinity matrix, to encode the edge-to-edge affinity between two graphs in the
off-diagonal elements. In this way, the two-graph matching between
$\mathcal{G}_{s}$ and $\mathcal{G}_{t}$ can be formulated as an edge-
preserving, quadratic assignment programming (QAP) problem [11]:

$\displaystyle\underset{\mathbf{v}}{\operatorname{argmax}}\
\mathbf{v}^{\top}M\mathbf{v}$ (1) $\displaystyle\text{ s.t. }\sum_{j\in
V_{t}}\mathbf{v}_{i,j}=1\;\forall i\in V_{s},\sum_{i\in
V_{s}}\mathbf{v}_{i,j}\leq 1\;\forall j\in V_{t},$

where $\mathbf{v}\in\\{0,1\\}^{nm\times 1}$ and $M\in\mathbb{R}^{nm\times
nm}$. The goal of graph matching is to establish a correspondence between two
attributed graphs, which minimizes the sum of local and geometric costs of
assignment between the vertices of the two graphs.

###  2.2 Position-Aware Node Embedding

Pixels at neighboring positions in the image convey similar semantic
information, therefore we need to distinguish them to resolve matching
ambiguities. Here, we refer to the nodes used as reference positions named as
anchors and propose an effective strategy to construct an anchor set $P$
consisting of all nodes in the graph, denoted by $P=V$, serving as a stable
reference for all nodes. To combine information from the nodes and the
anchors, we design the information aggregation mechanism as is shown in Fig.3.

Fig. 3: The position-aware embedding first aggregates the message of each
target node $h_{v}$ and each anchor point $h_{u_{i}}$ in the anchor set via
the aggregation function $I(v,u)$ based on the relative position coefficients
$q_{(}v,u)$. The representation matrix is then further aggregated using the
learnable function AGG and is finally transformed in a non-linear way to
obtain the feature output $h_{v}$ for the target node in the current layer.

Considering that the amount of effective information transmission to the nodes
by the anchors with different relative distances is different, we first
compute a relative position coefficient $q_{v,u}$ between a pair of nodes $v$
and $u$ as:

$q_{v,u}=\dfrac{e^{-d_{v,u}}}{\sum_{i\in V}{e^{-d_{v,i}}}},$ (2)

where $d_{v,u}$ is the shortest path distance between the two nodes. In
practice, to effectively reduce the time complexity of calculating the
shortest path and also to reduce interference from those anchor points that
are distant from the node under consideration, we demand that the maximum
shortest path distance does not exceed the ceiling value $r$. Otherwise, the
value of the coefficient is set to infinity. We continue to compute the
information aggregation function $I(v,u)^{(l)}$ for the $l$-th layer between
node $v$ and $u$ as:

$I(v,u)^{(l)}=q_{v,u}{\rm{CONCAT}}(h_{v}^{(l-1)},h_{u}^{(l-1)}),$ (3)

where $h_{v}^{(l-1)}$ and $h_{u}^{(l-1)}$ are the feature representation and
position representation in the $(l-1)$-th layer, and are combined through the
message aggregation function. We further aggregate the information from all
node-anchor pairs to obtaining a new representation in a high dimensional
space using the non-linear variation. This hidden representation can be
computed using the update

$h_{v}^{(l)}=\sigma({\rm{AGG}}(I(v,u)^{(l)}|\ \forall u\in V)W^{(l)}),$ (4)

where AGG is typically a permutation-invariant function (e.g., sum), and
$W^{(l)}$ is a learnable weight vector for the $l$-th layer.

###  2.3 Graph Feature Matching

Utilizing the proposed position-aware node embedding, the proposed scheme
encodes each node with position information into a high-level embedding space.
In this way, we can simplify the second-order affinity matrix of paired graphs
to a be a linear one with position-aware learning of the embedding. With the
final hidden representation of the source graph
$H_{s}\in\mathbb{R}^{n*F^{\prime}}$ and the target graph
$H_{t}\in\mathbb{R}^{m*F^{\prime}}$ to hand, we can obtain a soft
correspondence between these graphs with a node-wise affinity matrix through
the inner product of the embeddings of the two graphs being matched.
Specifically, to satisfy the condition that the original graph maps
injectively onto the target graph, we apply the Sinkhorn normalization [12] to
obtain rectangular doubly-stochastic correspondence matrices.

$S=\rm{sinkhorn}(H_{s}H_{t}^{T}).$ (5)

We further employ the cross-entropy loss function as the permutation loss
between the predicted permutation matrix and the ground truth:

$\displaystyle L=$
$\displaystyle-\sum_{i\in\mathcal{V}_{s},j\in\mathcal{V}_{t}}\left(\mathbf{S}_{i,j}^{gt}\log\mathbf{S}_{i,j}+\left(1-\mathbf{S}_{i,j}^{gt}\right)\log\left(1-\mathbf{S}_{i,j}\right)\right),$
(6)

where $\mathbf{S}^{gt}$ is the ground truth permutation matrix. Consequently,
the cross-entropy loss based on linear allocations can be learned end-to-end
no matter the number of nodes and edges in the graph.

Fig. 4:  Confusion matrix analysis of cross-category generalizations using the
Willow ObjectClass dataset. The $y$-axis represents categories of images used
for training the model, and the $x$-axis represents categories of samples used
to test. Among them, (a), (b), and (c) are the comparison of the
generalization ability of the aforementioned baseline in matching accuracy.
(d) represents the method proposed in this paper.

##  3 Experiments

###  3.1 Experimental Setting

Dataset. We use the Willow Object Class [13], consisting of 9963 images in
total and is used as a benchmark by many methods to measure the accuracy of
image classification and recognition algorithms. And IMC-PT-SparseGM [14]
datasets, containing 16 object categories and 25061 images, which gather from
16 tourist attractions around the world.
https://www.di.ens.fr/willow/research/graphlearning/.

Implementation Details. The experiments are conducted using two GeForce GTX
1080 Ti GPUs. We employ a batch size of 8 in training and evaluate the results
of our method after 2000 epochs for each iteration. We employ the Adam [15]
optimizer to train our models with a learning rate of $1\times 10^{-4}$. To
overcome the over-smoothing problem common to graph neural network models, we
adopt a two-layer graph embedding and restrict the degree of smoothing in our
experiments.

Metrics. We evaluate the graph matching capacity of the proposed method using
the matching accuracy, which is defined by
$accuracy=\sum{{\rm{AND}}(S_{i,j},{S_{i,j}^{gt}})/N}$ from [8, 16]. Note that
$S_{i,j}$ is the element of the predicted permutation matrix representing the
correspondence matching of node $i$ and node $j$ from two different graphs.
Similarly, $S_{i,j}^{gt}$ is the correspondence of the ground truth between
two nodes, and $N$ is the number of matching node pairs.

Table 1: Matching accuracy (%) on the Willow ObjectClass.

Method | Car | Duck | Face | M-bike | W-bottle | Mean
---|---|---|---|---|---|---
GMN[17] | 67.90 | 76.70 | 99.80 | 69.20 | 83.10 | 79.34
NHGM[18] | 86.50 | 72.20 | 99.90 | 79.30 | 89.40 | 85.50
CIE-H[19] | 82.20 | 81.20 | 100.00 | 90.00 | 97.60 | 90.20
PIA-GM[8] | 88.60 | 87.00 | 100.00 | 70.30 | 87.80 | 86.74
PCA-GM[8] | 87.60 | 83.60 | 100.00 | 77.60 | 88.40 | 87.44
IPCA-GM[16] | 90.40 | 88.60 | 100.00 | 83.00 | 88.30 | 90.06
PA-GM (ours) | 92.70 | 91.30 | 100.00 | 84.50 | 93.80 | 92.46

Table 2: Matching accuracy (%) on the IMC-PT-SparseGM.

Method | Reichstag | Sacre_coeur | St_peters_square | Mean
---|---|---|---|---
CIE-H[19] | 42.24 | 28.47 | 30.78 | 33.83
PIA-GM[8] | 71.46 | 41.31 | 42.64 | 51.80
PCA-GM[8] | 69.38 | 39.86 | 42.10 | 50.40
IPCA-GM[16] | 72.96 | 43.80 | 44.93 | 53.89
GANN-GM[20] | 76.02 | 44.15 | 50.49 | 56.89
PA-GM (ours) | 96.28 | 75.93 | 81.66 | 84.63

###  3.2 Graph Matching Results on Real-world Datasets

Table 1 and Table 2 report the overall comparison of the performance results
for graph matching accuracy. The CIE-H method excels on the M-bike and
W-bottle classes of the Willow ObjectClass dataset, while our method scores
highest in the two categories of accuracy, and the average accuracy. By
successfully incorporating the position information, our model achieves
excellent results in the average accuracy of Willow ObjectClass and IMC-PT-
SparseGM datasets. Based on these results, it can be concluded that our method
performs well on graph matching compared with existing methods.

###  3.3 Cross-category Generalization Study

To assess the robustness and generalization ability of our method on the
different object categories, we further conduct the cross-category
generalization study. Specifically, we train each of the five classes
separately, and then test the generalization ability of the validation model
for all five classes with the separate training classes. Comparing the
elements of Fig. 4, it is clear that the generalization ability of the model
framework that incorporates position information is superior to the
alternative methods. In addition, the matching accuracy for face data is
generally rather large. This may be related to the relatively simple
background of the data and less interference from noise.

##  4 Conclusion

In this paper, we propose a novel deep learning method for graph matching
based on node-wise embeddings between graphs that combine position-aware node
embeddings. During the experiments, the proposed method is compared with
alternative methods to demonstrate its robustness and effectiveness. Comparing
our methodology to existing methods on real-world datasets demonstrates its
state-of-the-art performance. Further improvements of graph matching will be
achieved with the use of different relative positions strategies in future
work.

## References

  * [1] Steven Gold and Anand Rangarajan,  “A graduated assignment algorithm for graph matching,”  IEEE Transactions on pattern analysis and machine intelligence, vol. 18, no. 4, pp. 377–388, 1996.
  * [2] Kexin Deng, Jie Tian, Jian Zheng, Xing Zhang, Xiaoqian Dai, and Min Xu,  “Retinal fundus image registration via vascular structure graph matching,”  International Journal of Biomedical Imaging, vol. 2010, 2010.
  * [3] Jiawei Zhang and S Yu Philip,  “Multiple anonymized social networks alignment,”  in 2015 IEEE International Conference on Data Mining. IEEE, 2015, pp. 599–608.
  * [4] Miao Wang, Yu-Kun Lai, Yuan Liang, Ralph R Martin, and Shi-Min Hu,  “Biggerpicture: data-driven image extrapolation using graph matching,”  ACM Transactions on Graphics, vol. 33, no. 6, 2014.
  * [5] Deepti Pachauri, Risi Kondor, and Vikas Singh,  “Solving the multi-way matching problem by permutation synchronization,”  in Advances in neural information processing systems. Citeseer, 2013, pp. 1860–1868.
  * [6] Junchi Yan, Shuang Yang, and Edwin R Hancock,  “Learning for graph matching and related combinatorial optimization problems,”  in Proceedings of the Twenty-Ninth International Joint Conference on Artificial Intelligence, IJCAI-20. International Joint Conferences on Artificial Intelligence Organization, 2020, pp. 4988–4996.
  * [7] Kun Xu, Liwei Wang, Mo Yu, Yansong Feng, Yan Song, Zhiguo Wang, and Dong Yu,  “Cross-lingual knowledge graph alignment via graph matching neural network,”  Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics, 2019.
  * [8] Runzhong Wang, Junchi Yan, and Xiaokang Yang,  “Learning combinatorial embedding networks for deep graph matching,”  in IEEE International Conference on Computer Vision, 2019, pp. 3056–3065.
  * [9] Karen Simonyan and Andrew Zisserman,  “Very deep convolutional networks for large-scale image recognition,”  3rd International Conference on Learning Representations, 2015.
  * [10] Hussam Qassim, Abhishek Verma, and David Feinzimer,  “Compressed residual-vgg16 cnn model for big data places image recognition,”  in 2018 IEEE 8th Annual Computing and Communication Workshop and Conference (CCWC). IEEE, 2018, pp. 169–175.
  * [11] Feng Zhou and Fernando De la Torre,  “Factorized graph matching,”  in 2012 IEEE Conference on Computer Vision and Pattern Recognition. IEEE, 2012, pp. 127–134.
  * [12] Richard Sinkhorn and Paul Knopp,  “Concerning nonnegative matrices and doubly stochastic matrices,”  Pacific Journal of Mathematics, vol. 21, no. 2, pp. 343–348, 1967\.
  * [13] Minsu Cho, Karteek Alahari, and Jean Ponce,  “Learning graphs to match,”  in International Conference on Computer Vision, 2013, pp. 25–32.
  * [14] Yuhe Jin, Dmytro Mishkin, Anastasiia Mishchuk, Jiri Matas, Pascal Fua, Kwang Moo Yi, and Eduard Trulls,  “Image matching across wide baselines: From paper to practice,”  International Journal of Computer Vision, pp. 517–547, 2021.
  * [15] Diederik P Kingma and Jimmy Ba,  “Adam: A method for stochastic optimization,”  International Conference on Learning Representations, 2015.
  * [16] Wang, Runzhong and Yan, Junchi and Yang, Xiaokang,  “Combinatorial learning of robust deep graph matching: an embedding based approach,”  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2020\.
  * [17] A. Zanfir and C. Sminchisescu,  “Deep learning of graph matching,”  in IEEE Conference on Computer Vision and Pattern Recognition, 2018, pp. 2684–2693.
  * [18] Runzhong Wang, Junchi Yan, and Xiaokang Yang,  “Neural graph matching network: Learning lawler’s quadratic assignment problem with extension to hypergraph and multiple-graph matching,”  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2021\.
  * [19] Tianshu Yu, Runzhong Wang, Junchi Yan, and Baoxin Li,  “Learning deep graph matching with channel-independent embedding and hungarian attention,”  in International Conference on Learning Representations, 2020.
  * [20] Runzhong Wang, Junchi Yan, and Xiaokang Yang,  “Graduated assignment for joint multi-graph matching and clustering with application to unsupervised graph matching network learning.,”  in NeurIPS, 2020.