# Generative Hypergraph Models and Spectral Embedding

Xue Gong School of Mathematics, University of Edinburgh, Edinburgh, EH9 3FD,
UK; The Maxwell Institute for Mathematical Sciences, Edinburgh, EH8 9BT, UK
<EMAIL_ADDRESS> Desmond J. Higham School of Mathematics, University
of Edinburgh, Edinburgh, EH9 3FD, UK<EMAIL_ADDRESS> Konstantinos
Zygalakis School of Mathematics, University of Edinburgh, Edinburgh, EH9 3FD,
UK<EMAIL_ADDRESS>
###### Abstract

Many complex systems involve interactions between more than two agents.
Hypergraphs capture these higher-order interactions through hyperedges that
may link more than two nodes. We consider the problem of embedding a
hypergraph into low-dimensional Euclidean space so that most interactions are
short-range. This embedding is relevant to many follow-on tasks, such as node
reordering, clustering, and visualization. We focus on two spectral embedding
algorithms customized to hypergraphs which recover linear and periodic
structures respectively. In the periodic case, nodes are positioned on the
unit circle. We show that the two spectral hypergraph embedding algorithms are
associated with a new class of generative hypergraph models. These models
generate hyperedges according to node positions in the embedded space and
encourage short-range connections. They allow us to quantify the relative
presence of periodic and linear structures in the data through maximum
likelihood. They also improve the interpretability of node embedding and
provide a metric for hyperedge prediction. We demonstrate the hypergraph
embedding and follow-on tasks—including quantifying relative strength of
structures, clustering and hyperedge prediction—on synthetic and real-world
hypergraphs. We find that the hypergraph approach can outperform clustering
algorithms that use only dyadic edges. We also compare several triadic edge
prediction methods on high school and primary school contact hypergraphs where
our algorithm improves upon benchmark methods when the amount of training data
is limited.

##  1 Introduction

A typical graph-based data set captures pairwise interactions between nodes.
There is growing interest in understanding higher-order, group-level,
interactions, with different paradigms being proposed [1, 2]. In this work, we
represent such interactions with a hypergraph formulation; here each hyperedge
involves two or more nodes. This framework is discussed in [3, 4, 5] and has
found application in real-world problems such as epidemic spread modelling
[6], image classification [7], and the study of biological networks [8].

A fundamental learning task on graph-based data is to embed nodes into a low-
dimensional Euclidean space [9, 10]. The learned embedding could be used in
follow-on tasks such as clustering, classification, and structure recovery.
There are various types of learning algorithms for graphs; some design and
analyze Laplacian matrices related to the graph [11], some solve maximum
likelihood problems associated with random graph models [12, 13], and others
involve more complex machine learning frameworks [14, 15].

In this work, we build on the use of spectral methods which derive node
embeddings from eigenvectors of a Laplacian matrix [16]. Such spectral
algorithms are popular, since they can be implemented efficiently on large
sparse graphs and they are backed up by accompanying consistency theory [17].
Two main approaches have also recently been proposed for spectral clustering
on hypergraphs. One approach is to employ higher-order Laplacian tensors [18].
Tensors in general contain richer information, however, their use can require
considerably more computational expense than matrix algorithms, and the
results can be difficult to visualize and interpret. A second approach is to
“flatten” the higher order information into a representative node-level
matrix. Some matrix-based approaches analyze the vertex-edge incidence matrix
associated with a random walk interpretation [9], other frameworks utilise
motif-based Laplacian matrices that could be generalized to various motifs and
time steps [10, 19]. The methodology that we develop here fits into this
second category by building a node-based matrix, using an intermediate step
that looks over all hyperedge dimensions in order to incorporate higher order
information.

A second aspect of our work is the connection between spectral methods and
random models. Many graph embedding [20], re-ordering [12], clustering [21],
and structure recovery [22, 23] techniques solve maximum-likelihood problems
on graphs assuming specific generative models. Besides their application in
these inverse problems, random graph models are useful inference tools for
quantifying structure, predicting new or missing links, and improving the
interpretability of learning algorithms by relating node embeddings to edge
probabilities [24, 25].

Many spectral algorithms are naturally related to optimization problems. This
is the case when the Laplacian matrix is Hermitian, so that its eigenvectors
are critical points of a quadratic form [26]. For example, spectral embedding
for undirected graphs using the standard combinatorial Laplacian is related to
minimizing the unnormalized cut [11, 27]. Furthermore, such optimization
formulations may lead to interesting random graph interpretations of spectral
algorithms. When the quadratic form can be expressed as the log-likelihood of
the graph under a suitable model, the optimization problem may be restated as
a node reordering. Such connections have been investigated for undirected
graphs [28] and directed graphs [25], and here we extend these ideas to the
hypergraph setting. In particular, we associate customized spectral embedding
algorithms with generative models that belong to a new class of range-
dependent random hypergraphs that encourages short-range connections between
nodes, generalizing existing graph models [12, 29]. These range-dependent
random hypergraphs offer flexibility that is not available in stochastic block
models [21] which require block sizes to be pre-specified or inferred.

The rest of the paper is structured as follows. Our notation is introduced in
Section 2. In Sections 3 and 4 we define the linear and periodic hypergraph
embedding algorithms and derive associated optimization problems. We propose
random models associated with the hypergraph embedding algorithms in Section
5, which leads to a model comparison workflow that quantifies the relative
strength of linear versus periodic structures. Numerical studies on synthetic
and real-world hypergraphs using the proposed models are presented in Section
6.

The main contributions of this work are as follows.

  * •

We propose new range-dependent generative models for hypergraphs that generate
linear and periodic cluster patterns.

  * •

We establish their connection with linear and periodic spectral embedding
algorithms.

  * •

We demonstrate on synthetic and real data that, after tuning model parameters
to the data, these models can quantify the relative strength of linear and
periodic structures.

  * •

We perform prediction of triadic hyperedges (triangles) using the proposed
linear model and show that it outperforms the existing average-score based
method [30] on synthetic hypergraphs, and also on high school and primary
school contact data when the amount of training data is limited.

##  2 Notation

We consider undirected, unweighted hypergraphs $G=(V,E)$ on the vertex set $V$
containing $n$ nodes and the hyperedge set $E$. We let $R\in\mathcal{R}$ be an
unordered set of nodes, where $\mathcal{R}$ denotes the collection of all such
sets. We use $\lvert R\rvert$ to denote the number of nodes in tuple $R$, that
is, its cardinality, and we assume $2\leq\lvert R\rvert\leq T$ for all $R\in
E$.

Let $A_{R}$ indicate the presence of a hyperedge, so that $A_{R}=1$ if $R\in
E$ and $A_{R}=0$ otherwise. We define the $t$-th order $n$ by $n$ adjacency
matrix $W^{[t]}$ such that $W^{[t]}_{ij}$ counts the number of hyperedges with
cardinality $t$ that contain distinctive nodes $i$ and $j$; hence,
$W^{[t]}_{ij}=\sum_{R\in\mathcal{R}:|R|=t}A_{R}\operatorname{\mathbbm{1}}{(i\in
R)}\operatorname{\mathbbm{1}}{(j\in R)}$ if $i\neq j$, and $W^{[t]}_{ij}=0$
otherwise, where $\operatorname{\mathbbm{1}}$ is the indicator function.
Similarly, we define the corresponding $t$-th order diagonal degree matrix
$D^{[t]}$ such that $D^{[t]}_{ii}=\sum_{j\in V}W^{[t]}_{ij}$, and the $t$-th
order Laplacian matrix $L^{[t]}=D^{[t]}-W^{[t]}$. We use
$\boldsymbol{\mathrm{i}}$ to denote $\sqrt{-1}$ and $1$ to denote the vector
in $\operatorname{\mathbb{R}}^{n}$ with all entries equal to one. We let
$\mbox{\boldmath$a$}^{\prime}$ represent the transpose of a real-valued vector
$a$ and let $\mbox{\boldmath$b$}^{H}$ denote the conjugate transpose of a
complex-valued vector $b$. We use $\mathcal{P}$ to denote the set of all
permutation vectors, that is, all vectors in $\operatorname{\mathbb{R}}^{n}$
that contain each of the integers $1,2,\ldots,n$. We will focus on one-
dimensional embedding. We let $x_{i}\in\operatorname{\mathbb{R}}$ be the
location to which node $i$ is assigned, and
$\mbox{\boldmath$x$}=[x_{1},x_{2},...,x_{n}]^{\prime}\in\operatorname{\mathbb{R}}^{n}$
.

##  3 Linear hypergraph embedding

Given a hypergraph, suppose we wish to find node embeddings
$\mbox{\boldmath$x$}\in\operatorname{\mathbb{R}}^{n}$ such that hyperedges
tend to contain nodes that are a small distance apart. To formalize this idea,
we can define a linear incoherence function
$I_{\text{lin}}(\mbox{\boldmath$x$},R)$ that sums up the squared Euclidean
distance between all nodes pairs in tuple $R$:

$I_{\text{lin}}(\mbox{\boldmath$x$},R)=\sum_{i,j\in R}(x_{i}-x_{j})^{2}.$ (1)

We may then define the total linear incoherence of the hypergraph,
$\eta_{\text{lin}}(G,\mbox{\boldmath$x$})$, by aggregating the linear
incoherence over all node tuples. Furthermore, we may wish to tune the weights
of hyperedges of different cardinalities through a coefficient $c_{\lvert
R\rvert}\geq 0$ for node tuple $R$; that is,

$\eta_{\text{lin}}(G,\mbox{\boldmath$x$})=\sum_{R\in\mathcal{R}}c_{\lvert
R\rvert}A_{R}I_{\text{lin}}(\mbox{\boldmath$x$},R).$ (2)

One justification for these tuning parameters $c_{\lvert R\rvert}$ is that
they allow us to avoid the case where high-cardinality hyperedges dominate the
expression. For example, we could choose $c_{t}=\frac{1}{t(t-1)}$ to balance
the contributions from hyperedges with different sizes. A suitable choice of
$c_{t}$ may also depend on the relative importance of hyperedges in the
application.

In Proposition 3.1 we show that the total linear incoherence maybe be written
as a quadratic form involving the hypergraph Laplacian matrix

$L=\sum_{t=2}^{T}c_{t}L^{[t]}.$ (3)

######  Proposition 3.1.

For any $\mbox{\boldmath$x$}\in\operatorname{\mathbb{R}}^{n}$ with $L$ defined
in (3), and $\eta_{\text{lin}}(G,\mbox{\boldmath$x$})$ defined in (2) we have

$\mbox{\boldmath$x$}^{\prime}L\mbox{\boldmath$x$}={\textstyle{{\frac{1}{2}}}}\eta_{\text{lin}}(G,\mbox{\boldmath$x$}).$
(4)

###### Proof.

It is straightforward to show that
$\mbox{\boldmath$x$}^{\prime}L^{[t]}\mbox{\boldmath$x$}=\sum_{i,j\in
V}x_{i}(D^{[t]}_{ij}-W^{[t]}_{ij})x_{j}={\textstyle{{\frac{1}{2}}}}\sum_{i,j\in
V}W^{[t]}_{ij}(x_{i}-x_{j})^{2}$. Therefore,

$\displaystyle\mbox{\boldmath$x$}^{\prime}L\mbox{\boldmath$x$}$
$\displaystyle=\sum_{t=2}^{T}\mbox{\boldmath$x$}^{\prime}c_{t}L^{[t]}\mbox{\boldmath$x$}$
$\displaystyle={\textstyle{{\frac{1}{2}}}}\sum_{t=2}^{T}c_{t}\sum_{i,j\in
V}W^{[t]}_{ij}(x_{i}-x_{j})^{2}$
$\displaystyle={\textstyle{{\frac{1}{2}}}}\sum_{t=2}^{T}c_{t}\sum_{i,j\in
V}\sum_{R\in\mathcal{R}:|R|=t}A_{R}\operatorname{\mathbbm{1}}{(i\in
R)}\operatorname{\mathbbm{1}}{(j\in R)}(x_{i}-x_{j})^{2}$
$\displaystyle={\textstyle{{\frac{1}{2}}}}\sum_{R\in\mathcal{R}}c_{\lvert
R\rvert}A_{R}\sum_{i,j\in
R}(x_{i}-x_{j})^{2}={\textstyle{{\frac{1}{2}}}}\eta_{\text{lin}}(G,\mbox{\boldmath$x$}).$

∎

We note that each Laplacian $L^{[t]}$ is symmetric and positive semi-definite
with smallest eigenvalue $0$.

######  Assumption 3.1.

We assume that the unweighted, undirected graph associated with the binarized
version of $L$ is connected. It then follows that $L$ has a single eigenvalue
equal to $0$ with all other eigenvalues positive. We further assume that there
is a unique smallest positive eigenvalue, $\lambda_{2}$. (The eigenvector
$\mbox{\boldmath$v$}^{[2]}$ corresponding to $\lambda_{2}$ is a generalization
of the classic Fiedler vector.)

In minimizing the total linear incoherence (2) we must avoid the trivial cases
where (a) all nodes are located arbitrarily close to the origin and (b) all
nodes are assigned to the same location. Hence it is natural to impose the
constraints $\|\mbox{\boldmath$x$}\|_{2}=1$ and
$\mbox{\boldmath$x$}^{\prime}\mbox{\boldmath$1$}=1$. It then follows from the
Rayleigh-Ritz Theorem [31, Theorem 4.2.2] that the quadratic form in
Proposition 3.1 is solved by $\mbox{\boldmath$x$}=\mbox{\boldmath$v$}^{[2]}$.
This leads us to Algorithm 1 below, which could also be considered as a
special case of the algorithm in [18] where the motifs considered are
hyperedges.

Result: Node embedding $\mbox{\boldmath$x$}\in\operatorname{\mathbb{R}}^{n}$

Input hyperedge adjacency matrices $W^{[2]},W^{[3]},...,W^{[T]}$;

Construct diagonal degree matrices $D^{[t]}_{ii}=\sum_{j\in V}W^{[t]}_{ij}$;

Construct $t$-th order Laplacians $L^{[t]}=D^{[t]}-W^{[t]}$;

Construct hypergraph Laplacian $L=\sum_{t=2}^{T}c_{t}L^{[t]}$;

Compute second smallest eigenvalue $\lambda_{2}$ and corresponding eigenvector
$\mbox{\boldmath$v$}^{[2]}$;

Embed nodes using $\mbox{\boldmath$x$}=\mbox{\boldmath$v$}^{[2]}$

Algorithm 1 Linear Hypergraph Embedding Algorithm

######  Remark 3.1.

Algorithm 1 could be extended to higher dimensional embeddings where node $i$
is assigned to $\mbox{\boldmath$x$}^{[i]}\in\operatorname{\mathbb{R}}^{d}$ for
$d>1$. In this case we could generalize (1) to

$I_{\text{lin}}(\mbox{\boldmath$x$},R)=\sum_{i,j\in
R}\lVert\mbox{\boldmath$x$}^{[i]}-\mbox{\boldmath$x$}^{[j]}\rVert_{2}^{2}.$
(5)

If we require the coordinate directions to be orthogonal, then the embedding
is found via the eigenvectors corresponding to the $d$-smallest non-zero
eigenvalues; see [27] for details in the graph case.

##  4 Periodic hypergraph embedding

In this section, we look at the periodic analogue of linear hypergraph
embedding. Here nodes are embedded into the unit circle rather than along the
real line. Such a periodic structure formed the basis of the classic “small
world” model of Watts and Strogatz [32]. Results in [24] showed that certain
real networks are better represented via this type of “wrap-around” notion of
distance. Hence, it is of interest to develop concepts that apply to the
hypergraph case.

We may position nodes on the unit circle by mapping them to phase angles
$\mbox{\boldmath$\theta$}=\\{\theta_{i}\\}_{i=1}^{n}\in[0,2\pi)$. We may then
use a periodic incoherence function to quantify the distance between node
pairs in the tuple $R$:

$I_{\text{per}}(\mbox{\boldmath$\theta$},R)=\sum_{i,j\in
R}|e^{\boldsymbol{\mathrm{i}}\theta_{i}}-e^{\boldsymbol{\mathrm{i}}\theta_{j}}|^{2}.$
(6)

Then the total periodic incoherence of the hypergraph becomes

$\eta_{\text{per}}(G,\mbox{\boldmath$\theta$})=\sum_{R\in\mathcal{R}}c_{\lvert
R\rvert}A_{R}I_{\text{per}}(\mbox{\boldmath$\theta$},R).$ (7)

In Proposition 4.1 below, we relate the total periodic incoherence to a
quadratic form involving the hypergraph Laplacian matrix (3).

######  Proposition 4.1.

Let $\mbox{\boldmath$\psi$}\in\operatorname{\mathbb{C}}^{n}$ be such that
$\psi_{j}=e^{\boldsymbol{\mathrm{i}}\theta_{j}}$. Then

$\mbox{\boldmath$\psi$}^{H}L\mbox{\boldmath$\psi$}={\textstyle{{\frac{1}{2}}}}\eta_{\text{per}}(G,\mbox{\boldmath$\theta$}).$
(8)

###### Proof.

We have

$\displaystyle\mbox{\boldmath$\psi$}^{H}L^{[t]}\mbox{\boldmath$\psi$}$
$\displaystyle=\mbox{\boldmath$\psi$}^{H}D^{[t]}\mbox{\boldmath$\psi$}-\mbox{\boldmath$\psi$}^{H}W^{[t]}\mbox{\boldmath$\psi$}$
$\displaystyle=\sum_{i\in V}e^{-\boldsymbol{\mathrm{i}}\theta_{i}}(\sum_{j\in
V}W^{[t]}_{ij})e^{\boldsymbol{\mathrm{i}}\theta_{i}}-\sum_{i,j\in
V}e^{-\boldsymbol{\mathrm{i}}\theta_{i}}W^{[t]}_{ij}e^{\boldsymbol{\mathrm{i}}\theta_{j}}$
$\displaystyle=\sum_{i,j\in
V}W^{[t]}_{ij}(1-e^{-\boldsymbol{\mathrm{i}}\theta_{i}}e^{\boldsymbol{\mathrm{i}}\theta_{j}})$
$\displaystyle={\textstyle{{\frac{1}{2}}}}\sum_{i,j\in
V}W^{[t]}_{ij}(2-e^{-\boldsymbol{\mathrm{i}}\theta_{i}}e^{\boldsymbol{\mathrm{i}}\theta_{j}}-e^{-\boldsymbol{\mathrm{i}}\theta_{j}}e^{\boldsymbol{\mathrm{i}}\theta_{i}})$
$\displaystyle={\textstyle{{\frac{1}{2}}}}\sum_{i,j\in
V}W^{[t]}_{ij}(e^{\boldsymbol{\mathrm{i}}\theta_{i}}-e^{\boldsymbol{\mathrm{i}}\theta_{j}})(e^{-\boldsymbol{\mathrm{i}}\theta_{i}}-e^{-\boldsymbol{\mathrm{i}}\theta_{j}})$
$\displaystyle={\textstyle{{\frac{1}{2}}}}\sum_{i,j\in
V}W^{[t]}_{ij}|e^{\boldsymbol{\mathrm{i}}\theta_{i}}-e^{\boldsymbol{\mathrm{i}}\theta_{j}}|^{2}.$

Then the proof may be completed in a similar way to the proof of Proposition
3.1. ∎

Appealing again to the Rayleigh–Ritz theorem [31, Theorem 4.2.2], the
quadratic form in (8) is minimized over all
$\mbox{\boldmath$\psi$}\in\operatorname{\mathbb{C}}^{n}$ with
$\|\mbox{\boldmath$\psi$}\|_{2}=1$ and
$\mbox{\boldmath$\psi$}^{H}\mbox{\boldmath$1$}=1$ by taking
$\mbox{\boldmath$\psi$}=\mbox{\boldmath$v$}^{[2]}$. However, this real-valued
eigenvector cannot be proportional to a vector with components of the form
$e^{\boldsymbol{\mathrm{i}}\theta_{j}}$. Hence, following the approach in [24]
we will use the heuristic of setting

$\theta_{i}=\text{angle}(v^{[2]}_{i}+\boldsymbol{\mathrm{i}}v^{[3]}_{i})\in[-\pi,\pi],$
(9)

defined as $v^{[2]}_{i}+\boldsymbol{\mathrm{i}}v^{[3]}_{i}=\lvert
v^{[2]}_{i}+\boldsymbol{\mathrm{i}}v^{[3]}_{i}\rvert\cdot e^{i\theta_{i}}$,
where $\mbox{\boldmath$v$}^{[3]}$ is an eigenvector corresponding to the next-
smallest eigenvalue of $L$. Such a heuristic assumption converts two real
eigenvectors into a complex vector, which gives an approximate solution to the
minimization problem. The resulting workflow is summarized in Algorithm 2.

We also note that for simple unweighted, undirected graphs, finding $\theta$
that minimizes $\eta_{\text{per}}(G,\mbox{\boldmath$\theta$})$ is equivalent
to the formulation proposed in [24]. This may be shown by letting
$u_{i}=\cos\theta_{i}$ and $z_{i}=\sin\theta_{i}$ and expanding (7) as

$\sum_{i\in V}\sum_{j\in
V}W^{[2]}_{ij}\left((u_{i}-u_{j})^{2}+(z_{i}-z_{j})^{2}\right),$

which simplifies to

$2\mbox{\boldmath$u$}^{T}(D^{[2]}-W^{[2]})\mbox{\boldmath$u$}+2\mbox{\boldmath$z$}^{T}(D^{[2]}-W^{[2]})\mbox{\boldmath$z$}.$
(10)

This is essentially equation (3.1) in [24], derived from a slightly different
perspective.

We then arrive at Algorithm 2 below.

Result: Node embedding
$\mbox{\boldmath$\theta$}=\\{\theta_{i}\\}_{i=1}^{n}\in[0,2\pi)$

Input adjacency matrices $W^{[2]},W^{[3]},...,W^{[T]}$;

Construct diagonal degree matrices $D^{[t]}_{ii}=\sum_{j\in V}W^{[t]}_{ij}$;

Construct $t$-th order Laplacians $L^{[t]}=D^{[t]}-W^{[t]}$;

Construct hypergraph Laplacian $L=\sum_{t=2}^{T}c_{t}L^{[t]}$;

Compute second and third smallest eigenvalues $\lambda_{2}$ and $\lambda_{3}$
and corresponding eigenvectors $\mbox{\boldmath$v$}^{[2]}$ and
$\mbox{\boldmath$v$}^{[3]}$;

Calculate phase angles
$\theta_{i}=\text{angle}(v^{[2]}_{i}+\boldsymbol{\mathrm{i}}v^{[3]}_{i})$;

Embed nodes using $\theta$

Algorithm 2 Periodic Hypergraph Embedding Algorithm

##  5 Generative hypergraph models

We now discuss a connection between the minimization of total incoherence and
generative models. Let us consider finding a node embedding
$\mbox{\boldmath$x$}\in\operatorname{\mathbb{R}}^{n}$ that minimizes a generic
total graph incoherence expression

$\eta(G)=\sum_{R\in\mathcal{R}}c_{\lvert
R\rvert}A_{R}I(\mbox{\boldmath$x$},R),$ (11)

for a non-negative incoherence function $I(\mbox{\boldmath$x$},R)$. We
consider the case where the $x_{i}\in\operatorname{\mathbb{R}}$ must take
distinct values from a discrete set $\\{\nu_{i}\\}_{i=1}^{n}$, where
$\nu_{i}\in\operatorname{\mathbb{R}}$; that is, we must have
$x_{i}=\nu_{p_{i}}$, where $p\in{\mathcal{P}}$ is a permutation vector. In the
linear case, this set may be the integers from $1$ to $n$ and in the periodic
case this set may be equally spaced angles in $[0,2\pi)$.

Now consider a random hypergraph model where each hyperedge involving node
tuple $R\in\mathcal{R}$ is generated independently with probability

$\textbf{P}(A_{R}=1)=f_{R}(\mbox{\boldmath$x$},R),$ (12)

for a function $f_{R}$ that takes values between 0 and 1. We have the
following connection.

######  Theorem 5.1.

Suppose $\mbox{\boldmath$x$}\in\operatorname{\mathbb{R}}^{n}$ is constrained
to take values from a discrete set such that $x_{i}=\nu_{p_{i}}$, where
$p\in{\mathcal{P}}$ is a permutation vector. Then minimizing the total
incoherence (11) over all such $x$ is equivalent to maximizing over all such
$x$ the likelihood that the hypergraph is generated by a model of the form
(12), where

$f_{R}(\mbox{\boldmath$x$},R)=\frac{1}{1+e^{\gamma c_{\lvert
R\rvert}I(\mbox{\boldmath$x$},R)}}$ (13)

for any positive $\gamma$.

###### Proof.

Using (12), the likelihood of the whole hypergraph is

$\displaystyle L(G)$
$\displaystyle=\prod_{R\in\mathcal{R}:A_{R}=1}f_{R}(\mbox{\boldmath$x$},R)\prod_{R\in\mathcal{R}:A_{R}=0}\left(1-f_{R}(\mbox{\boldmath$x$},R)\right)$
$\displaystyle=\prod_{R\in\mathcal{R}:A_{R}=1}\frac{f_{R}(\mbox{\boldmath$x$},R)}{1-f_{R}(\mbox{\boldmath$x$},R)}\prod_{R\in\mathcal{R}}\left(1-f_{R}(\mbox{\boldmath$x$},R)\right),$

which leads to the log-likelihood

$\displaystyle\ln(L(G))$
$\displaystyle=\sum_{R\in\mathcal{R}:A_{R}=1}\ln\left(\frac{f_{R}(\mbox{\boldmath$x$},R)}{1-f_{R}(\mbox{\boldmath$x$},R)}\right)+\sum_{R\in\mathcal{R}}\ln\left(\left(1-f_{R}(\mbox{\boldmath$x$},R)\right)\right).$
(14)

The second term on the right-hand side, which is the probability of the null
hypergraph, is independent of the the permutation. Hence, with (13),
maximizing the log-likelihood of the hypergraph is equivalent to minimizing

$\displaystyle\sum_{R\in\mathcal{R}:A_{R}=1}\ln\left(\frac{1-f_{R}(\mbox{\boldmath$x$},R)}{f_{R}(\mbox{\boldmath$x$},R)}\right)=\sum_{R\in\mathcal{R}}c_{\lvert
R\rvert}A_{R}\gamma I(\mbox{\boldmath$x$},R)=\gamma\,\eta(G).$ (15)

∎

######  Remark 5.1.

Theorem 5.1 could be extended to the case where node $i$ is assigned to
$\mbox{\boldmath$x$}^{[i]}\in\operatorname{\mathbb{R}}^{d}$ for $d>1$, and a
higher-dimensional incoherence function in (5) is considered. In this
scenario, we constrain
$\mbox{\boldmath$x$}^{[i]}\in\operatorname{\mathbb{R}}^{d}$ to take values
from a discrete set $\\{\mbox{\boldmath$\nu$}^{[i]}\\}_{i=1}^{n}$ where
$\mbox{\boldmath$\nu$}^{[i]}\in\operatorname{\mathbb{R}}^{d}$, such that
$\mbox{\boldmath$x$}^{[i]}=\mbox{\boldmath$\nu$}^{[p_{i}]}$ for a permutation
vector $p\in{\mathcal{P}}$. Then we could follow the same arguments as in
Theorem 5.1 to derive a model described by (12) and (13), where
$\mbox{\boldmath$x$}\in\operatorname{\mathbb{R}}^{n\times d}$ and
$I(\mbox{\boldmath$x$},R)=\sum_{i,j\in
R}\lVert\mbox{\boldmath$x$}^{[i]}-\mbox{\boldmath$x$}^{[j]}\rVert_{2}^{2}$.

For a hypergraph generated by model (13) the number of hyperedges connecting
the node tuple $R$ follows a Bernoulli distribution with probability
$1/1+e^{\gamma c_{\lvert R\rvert}I(\mbox{\boldmath$x$},R)}$. The log-odds of
the hyperedge decay linearly with the incoherence of the node tuple since

$\ln(f_{R}(\mbox{\boldmath$x$},R)/(1-f_{R}(\mbox{\boldmath$x$},R)))=-\gamma
c_{\lvert R\rvert}I(\mbox{\boldmath$x$},R),$

where the factor $\gamma c_{\lvert R\rvert}$ determines the decay rate. The
probability of a hyperedge is highest when all nodes overlap, i.e.,
$I(\mbox{\boldmath$x$},R)=0$, which gives a $1/2$ likelihood. If we generate
hyperedges in repeated trials for the node tuple $R$, the variance of the
number of hyperedges is $e^{c_{\lvert R\rvert}\gamma
I(\mbox{\boldmath$x$},R)}/(1+e^{c_{\lvert R\rvert}\gamma
I(\mbox{\boldmath$x$},R)})^{2}$. When $I(\mbox{\boldmath$x$},R)=0$, the
largest variance of $1/4$ is achieved. The expected total number of hyperedges
of the whole hypergraph $G$ can be expressed as

$\sum_{R\in\mathcal{R}}f_{R}(\mbox{\boldmath$x$},R)=\sum_{R\in\mathcal{R}}\frac{1}{1+e^{\gamma
c_{\lvert R\rvert}I(\mbox{\boldmath$x$},R)}}.$

We note that Theorem 5.1 introduces the extra scaling parameter $\gamma$. This
parameter plays no direct role in Algorithms 1 and 2. However, a value for
$\gamma$ is needed if we wish to compare the likelihoods of the two models
having inferred the embeddings. In principle, we may fit the parameter
$\gamma$ to a given hypergraph by matching the observed number of hyperedges
with their expectation. However, from a computational point of view, this is
rather challenging in general, since the computational complexity of the
expectation calculation is $\mathcal{O}(n^{T})$ when the maximum cardinality
of a considered hyperedge is $T$. Hence, given an embedding, in practice we
prefer to pick $\gamma$ by maximizing the likelihood, as described in the
following subsection.

### Model comparison

Under the assumption that a given hypergraph arose from a mechanism that
favours connections between “nearby” nodes (in some latent, unobservable
configuration), it is of interest to know whether a linear or periodic
distance provides a better description.We may address this question using a
model comparison approach. As in Sections 3 and 4, we consider one-dimensional
embeddings, such that both the linear and periodic version have $n+1$
parameters given the Laplacian coefficients: $n$ node embeddings and a decay
parameter $\gamma$. The node embeddings will be estimated from Algorithms 1
and 2. For any choice of $\gamma$, we may then calculate the corresponding
likelihood for each type of hypergraph, given the embedding. We may then
compare the models by reporting plots of likelihood versus $\gamma$ or by
reporting the maximum likelihood over all $\gamma$. We note that Theorem 5.1
states that node embeddings that minimize the incoherence also maximize the
graph likelihood under the given discrete constraints. We note that Algorithms
1 and 2 minimize linear and periodic incoherence after relaxing the discrete
constraints in Theorem 5.1 for computational feasibility. Such heuristics are
often used in discrete programming. Therefore instead of the exact maximum
likelihood, we get an estimated maximum likelihood. An overall workflow is
shown below in Algorithm 3.

Result: Comparison of possible graph structures

Input hypergraph $G$;

Compute linear embedding using Algorithm 1;

Compute periodic embedding using Algorithm 2;

Calculate maximum likelihood of linear model over $\gamma>0$ using (14), (13),
and (1) ;

Calculate maximum likelihood of periodic model over $\gamma>0$ using (14),
(13), and (6) ;

Compare likelihoods or report maxima

Algorithm 3 Model Comparison

##  6 Experiments

### Model Comparison

#### Synthetic Hypergraphs

In this section we test the performance on Algorithms 1, 2 and 3 in a
controlled setting. To do this, we generate hypergraphs with either linear or
periodic clustered structure using the proposed random model. For simplicity,
we only consider dyadic and triadic edges, although the experiments could be
extended to include higher-order hyperedges.

##### Linear hypergraph with clustered nodes

We first generate hypergraphs with $K$ planted clusters
$C_{1},C_{2},...,C_{K}$ of size $m$, and $n=mK$ nodes. We embed the nodes
using $x_{i}=\frac{2(l-1)}{K}+\sigma$ if $i\in C_{l}$, where
$\sigma\sim\mathrm{unif}(-a,a)$ is an additive uniform noise. Hyperedges are
then drawn randomly according to model (13) with the linear incoherence (1).

We note that, in practice, the embedding algorithms must choose values $c_{2}$
and $c_{3}$ in order to form the hypergraph Laplacian, and the model
comparison algorithm must choose a value for $\gamma$. We are therefore
interested in the sensitivity of the process with respect to $c_{2}$ and
$c_{3}$, and in the accuracy with which $\gamma$ can be estimated. We use
$c_{2}$, $c_{3}$ and $\gamma_{0}$ to denote parameters used by the generative
model to create the synthetic data; we also let $c^{*}_{2}$ and $c^{*}_{3}$
denote the corresponding parameters used in the spectral embedding algorithms
and let $\gamma^{*}$ represent an inferred value of $\gamma_{0}$. We choose
$c_{2}=1$ and $c_{3}=1/3$ so that the weight of a hyperedge is inversely
proportional to the number of node pairs involved. We let $m=50$, $K=5$,
$a=0.05$, and vary the decay parameter $\gamma_{0}$ from 0 to 10. Figure 1(a)
shows an example of the dyadic adjacency matrix, $W^{[2]}$, with
$\gamma_{0}=4$, where dots represent non-zeros. A corresponding triadic
adjacency matrix, $W^{[3]}$, is shown in Figure 1(b). In all our tests we
discard hypergraphs that do not satisfy Assumption 3.1.

For each synthetic hypergraph, we estimate the maximum log-likelihood assuming
a linear or a periodic structure using Algorithm 3. For each input decay
parameter $\gamma_{0}$, 40 hypergraphs are generated independently and the
average maximum log-likelihood is plotted in Figure 1(c). The shaded regions
represent the estimated $80\%$ confidence interval. In this case, the linear
model correctly achieves a higher average maximum log-likelihood. The tight
bound of the confidence interval suggests that the result is consistent across
random trials.

We then perform K-means clustering using the periodic and linear embeddings
assuming 5 clusters and plot the Adjusted Rand Index (ARI) [33, 34, 35] in
1(d). Here, a larger ARI indicates a better clustering result . The dotted
line shows the average over 40 independently trials for each $\gamma_{0}$
value and the shaded area is the estimated $80\%$ confidence interval. The
plot suggests that the clustering from the linear embedding outperforms the
clustering from the periodic embedding.

We are interested in the effect of parameters $c_{3}$ and $c_{3}^{*}$ that
control the weight of triadic edges in the random graph model and spectral
embedding algorithm respectively. To conduct an experiment, we fix the weight
of dyadic edges $c_{2}=1$, $c_{2}^{*}=1$, and decay parameter
$\gamma_{0}=\gamma^{*}=1$, while varying $c_{3}$ and $c_{3}^{*}$. The maximum
log-likelihood of the linear model (Figure 1(e)) and the ARIs using the linear
embedding (Figure 1(f)) are shown as heat-maps over $c_{3}$ and $c_{3}^{*}$.
Values are the average over 40 random trials. Overall, choosing
$c_{3}^{*}=c_{3}$, gives the highest maximum likelihood. Therefore, when the
true $c_{3}$ is not known, it could be estimated using a maximum likelihood
method. In terms of the clustering result we note that when $c_{3}$ is large,
for example, when $c_{3}>0.3$, using information from triadic edges by setting
$c_{3}^{*}>0$ achieves a better ARI than using only diadic edges, i.e.,
$c_{3}^{*}=0$. This is because a large $c_{3}$ encourages more triadic edges
to be formed within clusters, whereas a small $c_{3}$ leads to more triadic
edges between clusters. In general the larger the $c_{3}$, the less sensitive
the ARI is to the choice of $c_{3}^{*}$

(a)

(b)

(c)

(d)

(e) Linear log-likelihood of graph $\times 10^{-6}$

(f) Linear ARI

Figure 1: Model comparison experiments on synthetic linear hypergraphs. (a)
Black pixels in the dyadic adjacency matrix ($\gamma_{0}=4$) represent edges,
and they reveal 5 clusters in the diagonal blocks. (b) In the triadic
adjacency matrix, colors reflect the number of triangles shared between nodes
($\gamma_{0}=4$). (c) The linear embedding achieves higher log-likelihoods
than the periodic version for hypergraphs generated with different
$\gamma_{0}$. (d) Adjusted Rand Indices of K-means clustering based on the
linear and periodic embedding are plotted against $\gamma_{0}$. (e) Values in
the heatmap represent the maximum likelihoods of the linear model $\times
10^{-6}$ from Algorithm 3. The maxima are found along the diagonal when
$c_{3}^{*}=c_{3}$. (f) Values represent the Adjusted Rand Indices of K-means
clustering based on the linear embedding for different values of $c_{3}^{*}$
and $c_{3}$.

##### Periodic hypergraph with clustered nodes

To generate hypergraphs with periodic clusters, we use a node embedding based
on a vector of angles
$\boldsymbol{\theta}=(\theta_{1},\theta_{2},...,\theta_{n})^{T}$ in
$[0,2\pi)$, forming $K$ clusters $C_{1},C_{2},\ldots,C_{K}$ of size $m$. In
particular, we let $\boldsymbol{\theta}_{i}=\frac{2\pi(l-1)}{K}+\sigma$ if
$i\in C_{l}$ for $1\leq l\leq K$, where $\sigma\sim\mathrm{unif}(-a,a)$ is the
added noise. The hyperedges are generated using model (13), where the
incoherence function is defined in (6). We choose $a=0.05\pi$, $c_{2}=1$,
$c_{3}=1/3$ and vary the decay parameter $\gamma_{0}$. Examples of the dyadic
and triadic adjacency matrices with $\gamma_{0}=1$ are shown in Figure 2(a)
and 2(b).

Using the same approach as in the previous section, we compare the maximum
log-likelihood and ARIs assuming linear and periodic structures in Figure 2(c)
and 2(d). We see that the periodic model achieves a higher maximum, and on
average the periodic embedding produces higher ARIs.

Heat-maps in Figure 2(e) and 2(f) show results for different combinations of
$c_{3}$ and $c^{*}_{3}$ for the periodic embedding algorithm. These results
were generated in the same way as for Figures 1(c) and 1(d). Higher maximum
likelihoods are achieved near the diagonal where $c_{3}^{*}=c_{3}$, hence the
true parameters $c_{3}$ for the underlying hypergraph could be estimated using
the maximum likelihood method. As in the previous example, when $c_{3}\geq
0.3$, using the triadic edges ($c_{3}^{*}$) improves the ARI. When
$c_{3}<0.3$, increasing $c_{3}^{*}$ leads to an inferior clustering result.
However when $c_{3}\geq 0.3$, ARI becomes less sensitive to the choice of
$c_{3}^{*}$ as long as it is positive.

In summary, these tests indicate that the algorithms are able to correctly
distinguish between linear and periodic range-dependency when one such
structure is present in the data. We observed that setting $c_{3}^{*}>0$
improves the ARI when the triadic edges have a strong structural pattern; that
is, when $c_{3}$ is large. Moreover, when the true parameter $c_{3}$ is
unknown we recommend choosing $c_{3}^{*}$ based on a maximum likelihood
estimation, that is, finding the value $c_{3}^{*}$ that returns the largest
maxima in Algorithm 3. Such a choice also achieve reasonable ARIs in our
synthetic examples as shown in the diagonal entries in Figure 1(f) and 2(f).

(a)

(b)

(c)

(d)

(e) Log-likelihood of periodic model $\times 10^{-6}$

(f) Periodic ARI

Figure 2: Model comparison experiments on synthetic periodic hypergraphs. (a)
The bottom-left and top-right blocks in the dyadic adjacency matrix
($\gamma_{0}=1$) reveal a periodic clustered structure. (b) Similar periodic
clustered pattern is shown in the triadic adjacency matrix ($\gamma_{0}=1$).
(c) The periodic embedding achieves higher log-likelihoods than the linear
embedding for different values of $\gamma_{0}$. (d) Adjusted Rand Indices of
K-means clustering based on the linear and periodic embedding are plotted
against $\gamma_{0}$. (e) Values in the heatmap indicate the maximum
likelihoods of the periodic model $\times 10^{-6}$ from Algorithm 3 and the
maxima lie on the diagonal where $c_{3}^{*}=c_{3}$. (f) Values represent the
Adjusted Rand Indices of K-means clustering based on the periodic embedding
for different values of $c_{3}^{*}$ and $c_{3}$.

#### Real Hypergraphs

##### High School Contact Data

The high school contact data from [36] records the frequency of student
interaction. Students are represented as nodes, and contacts between two or
three students are registered as dyadic or triadic edges. We retrieved the
hypergraph from [21] containing 327 nodes, and only studied its dyadic and
triadic edges considering the computational complexity. We construct the
hypergraph Laplacian $L=c_{2}^{*}L^{[2]}+c_{3}^{*}L^{[3]}$ and perform linear
and periodic spectral embedding. For the linear embedding we map nodes into
3-dimensional Euclidean space using the eigenvectors corresponding to the
three smallest eigenvalues that are larger than $0.01$. We make this choice
because the eigenvector associated with the smallest non-zero eigenvalue has
only a few non-zero entries and leads to trivial clusters. We fix
$c_{2}^{*}=1$ and vary $c_{3}^{*}$ since only the relative weight
$c_{3}^{*}/c_{2}^{*}$ matters in node embedding.

The the maximum likelihoods and ARIs evaluated using various $c_{3}^{*}$ are
shown in the left column in Figure 3. The true clusters are defined by the
classes the students came from. Overall the periodic embedding achieves higher
likelihoods and ARIs despite the linear embedding involving more parameters.
Since linear clusters tend to have more marginalized groups that are far from
other clusters, our results may suggest a lack of marginalisation driven by
class membership.

We note that setting $c^{*}_{3}=0$ causes the algorithm to ignore triangles,
and hence to reduce to classical spectral clustering. For the linear
algorithm, we see that incorporating triadic edges by using a positive
$c^{*}_{3}$ can improve the ARI by up to around $0.09$. We note that in [21],
modularity maximization-based clustering achieved ARI=1 on the same data.
However, those methods have more parameters, which makes the ARI not directly
comparable.

##### Primary School Contact Data

The primary school contact hypergraph [21] is constructed from the contact
pattern between primary students from 10 classes [37]. Nodes represent
students or teachers, and hyperedges represent their physical contact. Each
node is labelled by the class of the students or as a teacher. The hypergraph
contains 242 nodes and 11 classes of labels. We extracted the dyadic and
triadic edges from the hypergraph and performed likelihood comparison and
clustering with four eigenvectors associated with the four smallest
eigenvalues that are greater than 0.01. The middle column in Figure 3 suggests
the periodic embedding achieves the maximum likelihood at $c_{3}^{*}=0.7$ and
overall performs better than the linear embedding in the clustering task.
These results may be related to the existence of a teacher group that connects
with all student groups. When we arrange dyadic and triadic adjacency matrices
by node classes, these connections will appear as off-diagonal entries. As we
have shown in Figure 2(a) and 2(b), the periodic model tends to produce more
off-diagonal connections than the linear model.

##### Senate Bills Data

In the senate bills hypergraph[38, 39, 21], nodes are US Congresspersons and
hyperedges are the sponsor and co-sponsors of Senate bills. There are in total
294 nodes, and each node is labelled as either Democrat or Republican. We
performed likelihood comparison and clustering with only the dyadic and
triadic edges. Since the node degree distribution is highly inhomogeneous, we
observe many trivial eigenvectors that are close to indicator functions. To
address this issue we trimmed off the top and bottom 2% nodes by node degree,
and use the eigenvector associated with the smallest eigenvalue that is
greater than 0.01. The linear and periodic models have similar maximum
likelihoods and clustering ARIs, as shown in the right column in Figure 3. In
contrast with previous examples, there are only two clusters present in this
data set. Hence the difference between the periodic and linear models, which
could be reflected in the connection (or disconnection) pattern between the
first and the last group, is less prominent.

High School Contact

Primary School Contact

Senate Bills

Figure 3: Model comparison experiments on the high school contact data (left),
the primary school contact data (middle), and the senate bill data (right).
Top row shows the maximum log-likelihoods of the linear and periodic embedding
for different values of $c_{3}^{*}$ where pentagrams indicate maxima. Bottom
row plots the Adjusted Rand Indices of K-means clustering based on the linear
and periodic embedding for various $c_{3}^{*}$.

### Hyperedge Prediction

Once the node embeddings are estimated from the spectral algorithms, the
probability of hyperedges may be computed from the proposed models. The
hyperedge probability can naturally serve as a score for hyperedge prediction.
We implement and test such triadic edge prediction on timestamped high school
contact data [30, 36] , primary school contact data [37, 30], and synthetic
linear hypergraphs. The results will be compared against approaches based on
average-scores proposed in [30]. Other hyperedge prediction methods include
feature-based prediction [40], model-based prediction[41], and machine
learning-based prediction [42].

For high school and primary school contact data, we used three and four
eigenvectors respectively corresponding to the smallest eigenvalues that are
greater than $0.01$ to be consistent with the previous section, and only
consider dyadic and triadic edges. The hyperedges are sorted by time stamps
and split into training and testing data. For example, an $80:20$
training/testing splitting ratio means we use the first $80\%$ of the
hyperedges to train the model and the last $20\%$ to test the predictions.
When the training ratio is low, the subgraph for training may be disconnected
thus violating Assumption 3.1. Therefore, we only consider nodes in the
largest connected component of the graph associated with the binarized version
of $L$ of the training subgraph, and test the prediction on the same set of
nodes. Note that in real data, the parameters $\gamma_{0}$, $c_{2}$ and
$c_{3}$ for the hypergraph model are unknown. We fix $c_{2}^{*}=1$ and choose
$c_{3}^{*}$ and $\gamma^{*}$ using a maximum likelihood estimate through a
grid search on the training data.

On the training set, we assign scores to each triplet using five methods: a
random score as baseline, hyperedge probability from the linear model,
arithmetic mean, harmonic mean, and geometric mean from [30]. On the test set,
we measure the prediction performance with the area under precision-recall
curve (AUC-PR) [43]. A Precision-Recall (PR) curve traces the Precision = True
Positive/ (True Positive + False Positive) and Recall = True Positive / (True
Positive + False Negative) for different thresholds. The AUC-PR is a measure
that balances both Precision and Recall where 1 means perfect prediction at
any threshold. Setting $c_{3}^{*}=0$ will assign probability of $0.5$ to all
triplets, which is equivalent to the random score approach if we break ties
randomly.

On the high-school contact data shown in Table 1, the harmonic and geometric
mean attain the highest AUC-PR for large amounts of training data, see, for
example the $80:20$ data split; while the linear model predictions achieve the
best results for small amounts of training data, as seen in the $20:80$ data
split. This could be because when training data is insufficient, there are
more unobserved ”missing” dyadic and triadic edges. In this case, the node
embedding algorithm can infer node proximity based on common neighbours. In
other words, it can place nodes with common neighbours nearby even if they
haven’t been directly linked before. On the other hand, the geometric and
harmonic mean will assign a score of zero to a triplet if none of the nodes
has been connected previously, and therefore will predict no triadic edges.

We also test the triadic edge prediction on the synthetic linear hypergraphs,
generated in the manner described in subsection 6, with $K=4$, $m=60$,
$\gamma=10$, $c_{2}=1$, and $c_{3}=0.3$, such that the clustered pattern
resembles the high-school contact data. We consider three eigenvectors
associated with the smallest eigenvalues that are greater than 0.01 for the
synthetic linear model. Since defining a periodic model with more than one
eigenvector is beyond the scope of this work, we only test the linear
hypergraphs. We randomly select a portion of the hyperedges as the training
set, while ensuring the sampled hypergraph is connected, and test the
performance on the rest of the hyperedges. The AUC-PR averaged over 20 random
hypergraphs is shown in Table 1. We observe that the linear model outperforms
random score and average scores for various training data sizes.

Data (Train:Test) | Random | Linear | Arith Mean | Geo Mean | Harm Mean
---|---|---|---|---|---
Highschool (80:20) | 5.3e-5 | 9.3e-4 | 9.4e-4 | 6.2e-3 | 6.0e-3
Highschool (60:40) | 1.2e-4 | 1.2e-3 | 2.1e-3 | 1.1e-2 | 1.1e-2
Highschool (20:80) | 2.6e-4 | 9.8e-3 | 3.9e-3 | 7.6e-3 | 7.4e-3
Primary School (80:20) | 3.2e-4 | 8.4e-3 | 3.3e-3 | 1.6e-2 | 1.7e-2
Primary School (60:40) | 1.1e-3 | 1.2e-2 | 9.9e-3 | 2.2e-2 | 2.2e-2
Primary School (20:80) | 1.6e-3 | 3.0e-2 | 1.6e-2 | 2.2e-2 | 2.1e-2
Linear Model (80:20) | 6.4e-3 | 1.6e-1 | 4.1e-2 | 7.8e-2 | 7.7e-2
Linear Model (60:40) | 1.3e-2 | 2.8e-1 | 8.5e-2 | 1.8e-1 | 1.8e-1
Linear Model (20:80) | 3.2e-2 | 3.0e-1 | 1.1e-1 | 1.4e-1 | 1.3e-1
Table 1: AUC-PR for triangle prediction on high-school contact data and
synthetic hypergraphs from the linear model. Highest values are indicated in
bold.

##  7 Conclusion

In this work we have developed new random models and embedding algorithms for
hypergraphs, and investigated their equivalence. In particular, we focused on
two spectral embedding algorithms customized for hypergraphs, which aim to
reveal linear and periodic structures, respectively. We also described random
hypergraph models associated with these algorithms, which allow us to quantify
the relative strength of linear and periodic structures based on maximum
likelihood. We demonstrated the model comparison approach on synthetic linear
and periodic hypergraphs, showing that the results are consistent with the
generating mechanism. When applied to high school and primary school contact
hypergraphs, the model comparison suggests the periodic structure is more
prominent. On this data set we also showed that the “spectral embedding plus
random hypergraph” approach gives a useful strategy for predicting new
hyperedges.

In future work, it would be interesting to investigate how these linear and
periodic hypergraph models compare with other versions that use alternative
assumptions, including those based on core-periphery [23] and stochastic block
model [13, 21] structures.

## Acknowledgements

X.G. acknowledges support of MAC-MIGS CDT under EPSRC grant EP/S023291/1.
D.J.H. was supported by EPSRC grant EP/P020720/1. K.C.Z. was supported by the
Leverhulme Trust grant 2020-310 and EPSRC grant EP/V006177/1.

## Data, code and materials

This research made use of public domain data that is available over the
internet, as indicated in the text. Code for the experiments is available at
https://github.com/OpalGX/hypergraph_model.

## Competing interests

The authors declare that they have no conflicts of interest.

## Authors’ contributions

X.G. carried out the numerical experiments and drafted the manuscript. All
authors contributed to the theoretical understanding, development of research,
and the writing of the manuscript. All authors gave final approval for
publication.

## References

  * [1] Bianconi G.  Higher order networks : An introduction to simplicial complexes.  Cambridge University Press; 2021.
  * [2] Lambiotte R, Rosvall M, Scholtes I.  From networks to optimal higher-order models of complex systems.  Nature Physics. 2019;15(4):313-20.
  * [3] Benson AR, Gleich DF, Leskovec J.  Higher-order organization of complex networks.  Science. 2016;353(6295):163-6.
  * [4] Benson AR, Gleich DF, Higham DJ.  Higher-order network analysis takes off, fueled by classical ideas and new data.  SIAM News. 2021.  Available from: https://sinews.siam.org/Details-Page/higher-order-network-analysis-takes-off-fueled-by-old-ideas-and-new-data.
  * [5] Torres L, Blevins AS, Bassett DS, Eliassi-Rad T.  The why, how, and when of representations for complex systems.  SIAM Review. 2021;63:435-85.
  * [6] Higham DJ, De Kergorlay HL.  Epidemics on hypergraphs: Spectral thresholds for extinction.  Proceedings of the Royal Society A. 2021;477(2252):20210232.
  * [7] Yu J, Tao D, Wang M.  Adaptive Hypergraph Learning and its Application in Image Classification.  IEEE Transactions on Image Processing. 2012;21(7):3262-72.
  * [8] Ramadan E, Tarafdar A, Pothen A.  A hypergraph model for the yeast protein complex network.  In: 18th International Parallel and Distributed Processing Symposium, 2004\. Proceedings.; 2004. p. 189.
  * [9] Schölkopf B, Platt J, Hofmann T.  Learning with Hypergraphs: Clustering, Classification, and Embedding.  In: Advances in Neural Information Processing Systems 19: Proceedings of the 2006 Conference; 2007. p. 1601-8.
  * [10] Rossi RA, Ahmed NK, Koh E, Kim S, Rao A, Abbasi-Yadkori Y.  A Structural Graph Representation Learning Framework.  In: Proceedings of the 13th International Conference on Web Search and Data Mining. WSDM 2020. New York, NY, USA: Association for Computing Machinery; 2020. p. 483–491.  Available from: https://doi.org/10.1145/3336191.3371843.
  * [11] Luxburg U.  A tutorial on spectral clustering.  Statistics and Computing. 2007;17(4):395-416.
  * [12] Grindrod P.  Range-dependent random graphs and their application to modeling large small-world Proteome datasets.  Physical Review E. 2002;66(6):066702/7-066702.
  * [13] Peixoto TP.  Ordered community detection in directed networks.  Phys Rev E. 2022 Aug;106:024305.  Available from: https://link.aps.org/doi/10.1103/PhysRevE.106.024305.
  * [14] Hamilton WL.  Graph Representation Learning.  Synthesis Lectures on Artificial Intelligence and Machine Learning. 2020;14(3):1-159.
  * [15] Hamilton WL, Ying R, Leskovec J.  Representation learning on graphs: Methods and applications.  IEEE Data Engineering Bulletin. 2017.
  * [16] Chung F.  Spectral Graph Theory.  Regional conference series in mathematics; no. 92. Providence, R.I.: American Mathematical Society; 1997.
  * [17] De Kergorlay HL, Higham DJ.  Consistency of anchor-based spectral clustering.  Information and Inference: A Journal of the IMA. 2022;11:801—822.
  * [18] Galuppi F, Mulas R, Venturello L.  Spectral theory of weighted hypergraphs via tensors.  Linear and Multilinear Algebra. 2022;0(0):1-31.  Available from: https://doi.org/10.1080/03081087.2022.2030659.
  * [19] Lucas M, Cencetti G, Battiston F.  Multiorder Laplacian for synchronization in higher-order networks.  Physical Review Research. 2020 Sep;2:033410.  Available from: https://link.aps.org/doi/10.1103/PhysRevResearch.2.033410.
  * [20] Hoff PD, Raftery AE, Handcock MS.  Latent space approaches to social network analysis.  Journal of the American Statistical Association. 2002;97(460):1090-8.
  * [21] Chodrow PS, Veldt N, Benson AR.  Generative hypergraph clustering: From blockmodels to modularity.  Science Advances. 2021;7(28):eabh1303.  Available from: https://www.science.org/doi/abs/10.1126/sciadv.abh1303.
  * [22] Cimini G, Mastrandrea R, Squartini T.  Reconstructing networks.  Cambridge University Press; 2021.
  * [23] Tudisco F, Higham DJ.  Core-periphery detection in hypergraphs.  arXiv preprint arXiv:220212769. 2022.
  * [24] Grindrod P, Higham DJ, Kalna G.  Periodic reordering.  IMA Journal of Numerical Analysis. 2010;30(1):195-207.
  * [25] Gong X, Higham DJ, Zygalakis K.  Directed network Laplacians and random graph models.  Royal Society Open Science. 2021;8(10):211144.  Available from: https://royalsocietypublishing.org/doi/abs/10.1098/rsos.211144.
  * [26] Lütkepohl H.  Handbook of Matrices.  Chichester: Wiley; 1996.
  * [27] Higham DJ, Kalna G, Kibble MJ.  Spectral clustering and its use in bioinformatics.  Journal of Computational and Applied Mathematics. 2007;204:25-37.
  * [28] Higham DJ.  Unravelling small world networks.  Journal of Computational and Applied Mathematics. 2003;158(1):61-74.
  * [29] Higham DJ.  Spectral reordering of a range-dependent weighted random graph.  IMA Journal of Numerical Analysis. 2005;25(3):443-57.
  * [30] Benson AR, Abebe R, Schaub MT, Jadbabaie A, Kleinberg J.  Simplicial closure and higher-order link prediction.  Proceedings of the National Academy of Sciences. 2018;115(48):E11221-30.
  * [31] Horn RA, Johnson CR.  Matrix Analysis.  Cambridge: Cambridge University Press; 1985.
  * [32] Watts DJ, Strogatz SH.  Collective dynamics of ‘small-world’ networks [Journal Article].  Nature. 1998;393(6684):440-2.  Available from: https://doi.org/10.1038/30918.
  * [33] Rand WM.  Objective criteria for the evaluation of clustering methods.  Journal of the American Statistical Association. 1971;66(336):846 850\.
  * [34] Hubert L, Arabie P.  Comparing partitions.  Journal of Classification. 1985;2(1):193-218.  Available from: https://doi.org/10.1007/BF01908075.
  * [35] McComb C.  Adjusted Rand Index.  GitHub. 2022.  Retrieved June 29, 2022.  Available from: https://github.com/cmccomb/rand_index.
  * [36] Mastrandrea R, Fournet J, Barrat A.  Contact patterns in a high school: a comparison between data collected using wearable sensors, contact diaries and friendship surveys.  PloS One. 2015;10(9):e0136497.
  * [37] Stehlé J, Voirin N, Barrat A, Cattuto C, Isella L, Pinton JF, et al.  High-Resolution Measurements of Face-to-Face Contact Patterns in a Primary School.  PLoS ONE. 2011;6(8):e23176.  Available from: https://doi.org/10.1371/journal.pone.0023176.
  * [38] Fowler JH.  Connecting the Congress: A Study of Cosponsorship Networks.  Political Analysis. 2006;14(04):456-87.  Available from: https://doi.org/10.1093/pan/mpl002.
  * [39] Fowler JH.  Legislative cosponsorship networks in the US House and Senate.  Social Networks. 2006 oct;28(4):454-65.  Available from: https://doi.org/10.1016/j.socnet.2005.11.003.
  * [40] Yoon S, Song H, Shin K, Yi Y.  How much and when do we need higher-order information in hypergraphs? a case study on hyperedge prediction.  In: Proceedings of The Web Conference 2020; 2020. p. 2627-33.
  * [41] Scholtes I.  When is a network a network? Multi-order graphical model selection in pathways and temporal networks.  In: Proceedings of the 23rd ACM SIGKDD international conference on knowledge discovery and data mining; 2017. p. 1037-46.
  * [42] Yadati N, Nitin V, Nimishakavi M, Yadav P, Louis A, Talukdar P.  NHP: Neural hypergraph link prediction.  In: Proceedings of the 29th ACM International Conference on Information & Knowledge Management; 2020. p. 1705-14.
  * [43] Davis J, Goadrich M.  The relationship between Precision-Recall and ROC curves.  In: Proceedings of the 23rd International Conference on Machine Learning; 2006. p. 233-40.


Transport via Partial Galois Connections and Equivalences (Extended Version)
Transport via Partial Galois Connections and Equivalences

Anon. AuthorKevin Kappelmann0000-0003-1421-6497

Anon. PlaceTechnical University of Munich, Boltzmannstrasse 3, Garching 85748, Germany,
anon<EMAIL_ADDRESS>
§ APPLICATION EXAMPLES

As all our results are formalised in Isabelle/HOL,
we can directly use them to manually transport terms in said environment.
But that would be rather tiresome.

We thus implemented a prototype in Isabelle/ML to
automate transports.

The Prototype

The method $\transporttermprover$
uses registered base equivalences,
along with the closure theorems from <ref>,
to construct more complex equivalences.
The prototype is currently
restricted to equivalences on partial equivalence relations (PERs) for pragmatic reasons:

their closure theorems have fewer
assumptions and are hence simpler to apply.

Providing automation for weaker equivalences is future work.
The current prototype also does not build
composition closures (<ref>)
and automates only a fragment of
dependent function relators for simplicity reasons.
Again, these extensions are future work.

The prototype provides a command $\transportterm$.
As input, it takes a term $t\holhasty\alphaty$ (the term to be transported)
and two optional target relations
This is unlike other transport frameworks [29, 26, 9, 13],
which only take the term $t\holhasty \alphaty$ and a target type $\betaty$.
This design decision is crucial since we
can neither assume a unique correspondence
between types and target relations in practice
(cf. <ref>),
nor can we express dependencies in types,
but we express them using dependent relators
(cf. <ref>).
The command then opens two goals.
The first one asks for an equivalence
the second one for a proof that $\indom\app \inflerel{L}\app t$.
On success, it registers a new term $t'$ and a theorem that $t\galrel{L} t'$.
It also registers a second theorem where the relator $\infgalrel{L}$ has been rewritten
to its desired form
as described in <ref>.

The following examples are best explored interactively in our formalisation.
We define the restricted equality relation on predicates as
$x =_P y\define P\app x\land x = y$
the restricted equality relation on sets as
$x =_S y\define x\in S\land x = y$.

It is easy to transport the list and set examples from <ref>.
We just have to prove the equivalence between
$\listfsetrelL\app xs\app xs'\define \listfsetrel\app xs\app (\tofset\app xs')$
$(=) \holhasty \relty{\fsetty{\natty}}{\fsetty{\natty}}$
and invoke our prototype on $\maxlist$:

The $\perintrotag$ tag

is used by $\transporttermprover$
to discharge the closure theorems' side conditions.

$\transportterm$ registers the theorem
$\parenths[\big]{\funrel{\listfsetrel}{(=)}}\app \maxlist\app \maxfset$
and the definition
$\maxfset\app s\define \maxlist\app (\tolistfset\app s)$
as a result.

We can also readily transport in the opposite direction
or use $\settyconst$s rather than $\fsettyconst$s
if we
define $\listsetrelL\app xs\app xs'\define \listsetrel\app xs\app (\toset\app xs')$:

' t

As noted in <ref>, transporting subtractions
$i_1 -_\intty i_2$ from $\intty$ to $\natty$
requires a dependency $i_1 \geq i_2$.
We model this dependency using dependent function relators.
We first define $\intspos\define (=_{(\le) 0})$ and then
proceed as usual:

(-_)t (-_)

L []i_1_i_2_i_1 ≥i_2

R []n_1_(=)n_2_(=)n_1 ≥n_2(=)
Similarly, operations on datatypes may only conditionally be transportable.
For example, we may only transport the index operator
$\inflistidx\holhasty \listty{\alphaty}\purefun\natty\purefun\alphaty$
to the type of immutable arrays ($\iarrayty{\alphaty}$)
if the index is not out of bounds.
In the following, let $S$ be an arbitrary partial equivalence relation:

L []xs_Si_(=)i < xsS

R []arr_Si_(=)i < arrS

Isabelle/Set [14]
is a

set-theoretic environment in Isabelle/HOL.
Its type of sets is called $\mysetty$.
Isabelle/Set provides a set-extension mechanism:
As input, it takes two sets $A\holhasty\mysetty$ and $B\holhasty\mysetty$
and an injection from $A$ to $B$.
It then creates a new set $B'\supseteq A$ together with a bijection
between $B$ and $B'$
with mutual inverses $l,r \holhasty\mysetty\purefun\mysetty$.
This mechanism is used to enforce subset relationships.
For instance, it first uses a construction of the integers $\intty\holhasty\mysetty$
where $\natty\not\subseteq\intty$.
It then uses the set-extension mechanism to create a copy $\intty'\supseteq \natty$ with inverses $l,r$.
Doing so necessitates a manual transport of all definitions from
$\intty$ to $\inttycopy$.
Using ,
it is possible to automate this process:

(+_)t (+_)
L [](=_)(=_)(=_)

R [](=_)(=_)(=_)

(-_)t (-_)
L [](=_)(=_)(=_)

R [](=_)(=_)(=_)

Note that all constants $(+_\intty),(+_\inttycopy),(-_\intty),(-_\inttycopy)$
are of the same type $\mysetty\purefun\mysetty\purefun\mysetty$.
This stresses the point that users
must be able to specify target relations and not just target types.

§ RELATED WORK

Transport in Proof Assistants
Our work was chiefly inspired by Isabelle's

Lifting package [13, 17],
which transports terms via partial quotient types.

All closure theorems in this work

generalise the ones in [13].

Besides this source of inspiration,

the theory of automated transports
has seen prolific work in recent years:

Tabareau et al. 28
proved a strengthened relational parametricity result,
called univalent parametricity,
for the Calculus of Inductive Constructions.

Their approach ensures
that all relations are compatible with type equivalences.
One can then use univalence [33] to seamlessly transport
terms between related types.
The framework is implemented using Coq's
typeclass mechanism [27].

Tabareau et al. 29 extended their
work to integrate what they call
“white-box transports”.
White-box transports
structurally rewrite a term $t$ to $t'$ using
user-specified correspondences.
In contrast,
“black-box transports”

transport $t$ without looking at its syntactic structure.

For instance,
given an equivalence between unary and binary numbers
black-box transporting the term $0 +_{\natty} 0$ results
in $l\app (0 +_{\natty} 0)$.
In contrast,
given correspondences
between the functions ${(+)}_\natty, {(+)}_{\constfont{Bin}}$
and constants $0, 0_{\constfont{Bin}}$,
white-box transporting the term results
in $0_{\constfont{Bin}} +_{\constfont{Bin}} 0_{\constfont{Bin}}$.

These modes can also be mixed:
given just the equivalence
and correspondence between ${(+)}_{\natty},{(+)}_{\constfont{Bin}}$,
we obtain $(l\app 0) +_{\constfont{Bin}} (l\app 0)$.
Isabelle's Lifting package also supports white-box transports
via the $\transfer$ method [17].

While our work is concerned with
black-box transports,
our prototype also contains
experimental support for white-box transports.
This integration will be further polished in future work.

Angiuli et al. 1
establish representation independence results
in Cubical Agda [32].
Their approach

applies to a restricted variant of
quasi-partial equivalence relations [16].
Essentially, they quotient two types by a given correspondence to obtain
a type equivalence between the quotiented types.

Dagand et al. 9, 8 introduce what they call
“type-theoretic partial Galois connections”,
which are essentially partial type equivalences
on an enriched
$\alphaty\app \typefont{option}$ type.
They allow for partiality on one
side of the equivalence but not the other.
Their framework is designed for effective
program extraction and implemented
using Coq's typeclass mechanism.

Ringer et al. 26
developed a Coq plugin

to transport proof terms via type equivalences
for inductive types.
Their theory shares similarities
with [29, 28],
but it directly transforms proof terms.

This way, one can remove all references to the old datatype
once the proof terms have been transported to the new target type.

This is not readily achievable using other mentioned frameworks,
including ours.

Type equivalences enjoy the property
of having total and mutually inverse transport functions.
This is not the case for partial Galois connections,
which makes the transport of proofs harder.
For example, the parametricity law for equality
$\parenths[\big]{\funrel{T}{\funrel{T}{(\holiff)}}}\app (=)\app (=)$
holds only if $T$ is left-unique and injective.
This is the case if $T$ is described by a type equivalence
but not in general by a Galois connection.
Kunčar 17 provides parametricity rules
for all prominent logical connectives.
These rules also apply to our setting and
will be crucial when we polish the integration
of white-box transports in our prototype.

The works mentioned above
all transport terms via certain notions of equivalences.
But there are also other approaches,
particularly in the field of data refinement.
An example is the CoqEAL framework [4],
which automatically derives parametricity results using typeclass search.
Another one is Isabelle's Autoref framework [18],
which derives relational parametricity results using white-box transports.
The core inspiration in both cases goes back to
[21, 25, 34].
A comprehensive comparison of these frameworks
can be found in 19.

Galois Connections in Computer Science
Galois connections are fundamental

in the field of abstract interpretation.
Cousot and Cousot's recent book 5
provides an overview of their applications.

The closure of Galois connections under non-dependent function relators

goes back to at least [6].

We generalised this result
to partial Galois connections and dependent function relators
in <ref>.
Most work in abstract interpretation
does not consider partially defined Galois connections
and assumes partial orderings on relations.
The work of Miné 20
is an exception,
allowing for partiality on one side of the connection but not the other.
Darais and Van Horn 10 formalise Galois connections
constructively and apply it to tasks in abstract interpretation.
An early

application of Galois connections

was by Hartmanis and Stearns 12.

Though they did not use Galois connections,
they introduced an equivalent notion of pair algebras [11].
Our Galois relator

describes the pair algebra induced by a Galois connection.

§ CONCLUSION AND FUTURE WORK

We explored existing notions of equivalences used for automatic transport.
Based on this exploration,
we identified a set of minimal expectations when transporting terms via equivalences.
This essence led us to introduce a new class of equivalences,
namely partial Galois connections.
Partial Galois connections generalise (standard) Galois connections
and apply to relations
that are only defined on subsets of their types.
We derived closure conditions for partial Galois connections and equivalences,
and typical order properties
under (dependent) function relators, relators for (co)datatypes,
and composition.
Our framework applies to simple type theory and –
unlike prior solutions for simple type theory –
can handle inter-argument dependencies.
We implemented a prototype in Isabelle/HOL based on our results.
The prototype needs to be further polished,
but it can already handle relevant examples that
are out of scope for existing tools.

Future work

As our theory subsumes the one of Isabelle's Lifting package,
one goal is to replace the package by a more general tool.
To this end,
we have to integrate our results into Isabelle's (co)datatypes package [2],
extend our prototype to automate the construction of compositions,
and polish the support of white-box transports (cf. <ref>).

based on our formalisation insights,
we conjecture that one can adopt our theory
to constructive logics,
but only a formalisation in a constructive prover
will give a definite answer.

§.§.§ Acknowledgements

authors thankauthor thanks
the anonymous reviewers of this
and a previous submission
for their valuable feedback
and Mohammad Abdulaziz and Tobias Nipkow
for their comments on a draft of this paper.

[1]
Angiuli, C., Cavallo, E., Mörtberg, A., Zeuner, M.: Internalizing
Representation Independence with Univalence. Proc. ACM Program. Lang.
5(POPL) (jan 2021). 10.1145/3434293

[2]
Blanchette, J.C., Hölzl, J., Lochbihler, A., Panny, L., Popescu, A.,
Traytel, D.: Truly Modular (Co)datatypes for Isabelle/HOL. In: Klein, G.,
Gamboa, R. (eds.) Interactive Theorem Proving. pp. 93–110. Springer
International Publishing, Cham (2014). 10.1007/978-3-319-08970-6_7

[3]
Church, A.: A Formulation of the Simple Theory of Types. The Journal of
Symbolic Logic  5(2),  56–68 (1940). 10.2307/2266170

[4]
Cohen, C., Dénès, M., Mörtberg, A.: Refinements for Free! In:
Gonthier, G., Norrish, M. (eds.) Certified Programs and Proofs. pp. 147–162.
Springer International Publishing, Cham (2013).

[5]
Cousot, P.: Principles of Abstract Interpretation. MIT Press (2021)

[6]
Cousot, P., Cousot, R.: Static Determination of Dynamic Properties of
Recursive Procedures. In: Neuhold, E. (ed.) IFIP Conf. on Formal
Description of Programming Concepts, St-Andrews, N.B., CA. pp. 237–277.
North--Holland (1977)

[7]
Cousot, P., Cousot, R.: Abstract Interpretation Frameworks. Journal of Logic
and Computation  2(4),  511–547 (08 1992).

[8]
Dagand, P.E., Tabareau, N., Tanter, E.: Partial Type Equivalences for Verified
Dependent Interoperability. SIGPLAN Not.  51(9),  298–310 (sep
2016). 10.1145/3022670.2951933

[9]
Dagand, P.E., Tabareau, N., Tanter, E.: Foundations of Dependent
Interoperability. Journal of Functional Programming  28 (2018).

[10]
Darais, D., Van Horn, D.: Constructive Galois Connections. Journal of
Functional Programming  29 (2019). 10.1017/S0956796819000066

[11]
Derderian, J.C.: Galois Connections and Pair Algebras. Canadian Journal of
Mathematics  21,  498–501 (1969). 10.4153/CJM-1969-056-x

[12]
Hartmanis, J., Stearns, R.: Pair Algebra and Its Application to Automata
Theory. Information and Control  7(4),  485–507 (1964).

[13]
Huffman, B., Kunčar, O.: Lifting and Transfer: A Modular Design for
Quotients in Isabelle/HOL. In: Gonthier, G., Norrish, M. (eds.) Certified
Programs and Proofs - Third International Conference, CPP 2013, Melbourne,
VIC, Australia, December 11-13, 2013, Proceedings. Lecture Notes in Computer
Science, vol. 8307, pp. 131–146. Springer (2013).

[14]
Kappelmann, K., Josh, C., Krauss, A.: Isabelle/Set (2023),

[15]
Kappelmann, Kevin: Transport via Partial Galois Connections
and Equivalences. In: Hur, Chung-Kil (ed.) Asian Symposium on
Programming Languages and Systems. pp. 225–245. Springer, Singapore
(2023). 10.1007/978-981-99-8311-7_11

[16]
Krishnaswami, N.R., Dreyer, D.: Internalizing Relational Parametricity in the
Extensional Calculus of Constructions. In: Rocca, S.R.D. (ed.) Computer
Science Logic 2013 (CSL 2013). Leibniz International Proceedings in
Informatics (LIPIcs), vol. 23, pp. 432–451. Schloss
Dagstuhl–Leibniz-Zentrum fuer Informatik, Dagstuhl, Germany (2013).

[17]
Kunčar, O.: Types, Abstraction and Parametric Polymorphism in
Higher-Order Logic. Ph.D. thesis, Technische Universität München (2016)

[18]
Lammich, P.: Automatic Data Refinement. In: Blazy, S., Paulin-Mohring, C.,
Pichardie, D. (eds.) Interactive Theorem Proving. pp. 84–99. Springer Berlin
Heidelberg, Berlin, Heidelberg (2013). 10.1007/978-3-642-39634-2_9

[19]
Lammich, P., Lochbihler, A.: Automatic Refinement to Efficient Data
Structures: A Comparison of Two Approaches. Journal of Automated Reasoning
63(1),  53–94 (Jun 2019). 10.1007/s10817-018-9461-9

[20]
Miné, A.: Weakly Relational Numerical Abstract Domains. Theses, Ecole
Polytechnique X (Dec 2004),

[21]
Mitchell, J.C.: Representation Independence and Data Abstraction. In:
Proceedings of the 13th ACM SIGACT-SIGPLAN Symposium on Principles of
Programming Languages. p. 263–276. POPL '86, Association for Computing
Machinery, New York, NY, USA (1986). 10.1145/512644.512669

[22]
Moura, L.d., Kong, S., Avigad, J., van Doorn, F., von Raumer, J.: The Lean
Theorem Prover (System Description). In: Felty, A.P., Middeldorp, A. (eds.)
Automated Deduction - CADE-25. pp. 378–388. Springer International
Publishing, Cham (2015). 10.1007/978-3-319-21401-6_26

[23]
Moura, L.d., Ullrich, S.: The Lean 4 Theorem Prover and Programming Language.
In: Platzer, A., Sutcliffe, G. (eds.) Automated Deduction – CADE 28. pp.
625–635. Springer International Publishing, Cham (2021).

[24]
Nipkow, T., Wenzel, M., Paulson, L.C.: Isabelle/HOL: A Proof Assistant for
Higher-Order Logic. Springer-Verlag, Berlin, Heidelberg (2002).

[25]
Reynolds, J.C.: Types, Abstraction and Parametric Polymorphism. In: Mason,
R.E.A. (ed.) Information Processing 83, Proceedings of the IFIP 9th World
Computer Congress, Paris, France, September 19-23, 1983. pp. 513–523.
North-Holland/IFIP (1983)

[26]
Ringer, T., Porter, R., Yazdani, N., Leo, J., Grossman, D.: Proof Repair
across Type Equivalences. In: Proceedings of the 42nd ACM SIGPLAN
International Conference on Programming Language Design and Implementation.
p. 112–127. PLDI 2021, Association for Computing Machinery, New York, NY,
USA (2021). 10.1145/3453483.3454033

[27]
Sozeau, M., Oury, N.: First-Class Type Classes. In: Mohamed, O.A.,
Muñoz, C., Tahar, S. (eds.) Theorem Proving in Higher Order Logics. pp.
278–293. Springer Berlin Heidelberg, Berlin, Heidelberg (2008).

[28]
Tabareau, N., Tanter, E., Sozeau, M.: Equivalences for Free: Univalent
Parametricity for Effective Transport. Proc. ACM Program. Lang.
2(ICFP) (jul 2018). 10.1145/3236787

[29]
Tabareau, N., Tanter, E., Sozeau, M.: The Marriage of Univalence and
Parametricity. J. ACM  68(1) (jan 2021). 10.1145/3429979

[30]
Traytel, D., Popescu, A., Blanchette, J.C.: Foundational, Compositional
(Co)datatypes for Higher-Order Logic: Category Theory Applied to Theorem
Proving. In: 2012 27th Annual IEEE Symposium on Logic in Computer Science.
pp. 596–605 (2012). 10.1109/LICS.2012.75

[31]
Univalent Foundations Program, T.: Homotopy Type Theory: Univalent
Foundations of Mathematics. <https://homotopytypetheory.org/book>,
Institute for Advanced Study (2013)

[32]
Vezzosi, A., Mörtberg, A., Abel, A.: Cubical Agda: A Dependently Typed
Programming Language with Univalence and Higher Inductive Types. Proc. ACM
Program. Lang.  3(ICFP) (jul 2019). 10.1145/3341691

[33]
Voevodsky, V.: The equivalence axiom and univalent models of type theory
(2010). 10.48550/ARXIV.1402.5556

[34]
Wadler, P.: Theorems for free! In: Proceedings of the Fourth International
Conference on Functional Programming Languages and Computer Architecture. p.
347–359. FPCA '89, Association for Computing Machinery, New York, NY, USA
(1989). 10.1145/99370.99404


# Characteristic interaction potential of black hole molecules from the
microscopic interpretation of Ruppeiner geometry

Shao-Wen Wei1,2,3111Corresponding author. E-mail<EMAIL_ADDRESS>Yu-Xiao
Liu1,2,3222E-mail<EMAIL_ADDRESS>and Robert B. Mann4 333E-mail:
<EMAIL_ADDRESS>1Lanzhou Center for Theoretical Physics, Key Laboratory of
Theoretical Physics of Gansu Province, School of Physical Science and
Technology, Lanzhou University, Lanzhou 730000, People’s Republic of China,
2Institute of Theoretical Physics $\&$ Research Center of Gravitation, Lanzhou
University, Lanzhou 730000, People’s Republic of China,
3Academy of Plateau Science and Sustainability, Qinghai Normal University,
Xining 810016, P. R. China
4Department of Physics and Astronomy, University of Waterloo, Waterloo,
Ontario, Canada, N2L 3G1

###### Abstract

Ruppeiner geometry has been found to be a novel promising approach to uncover
the microstructure of fluid systems and black holes. In this work, combining
with the micro model of the Van der Waals fluid, we shall propose a first
microscopic interpretation for the empirical observation of Ruppeiner
geometry. Then employing the microscopic interpretation, we disclose the
potential microstructure for the anti-de Sitter black hole systems. Of
particular interest, we obtain the microscopic interaction potentials for the
underlying black hole molecules. This significantly strengthens the study
towards to the black hole nature from the viewpoint of the thermodynamics.

###### pacs:

04.70.Dy, 04.60.-m, 05.70.Ce

##  I Introduction

Classically a black hole is a purely absorbing object; it cannot emit
particles. However Hawking showed that when quantum mechanical effects are
considered near the horizon, black holes can indeed emit particles, with an
average energy flux given by Hawking

$\frac{dE}{dt}=\int_{0}^{\infty}\frac{\Gamma\omega}{e^{\omega\beta}-1}d\omega,$
(1)

where $\Gamma$ is the greybody factor of the radiation. This suggests that a
black hole emits particles at a temperature $T_{\rm H}=1/\beta$ known as the
Hawking temperature, which depends on the black hole mass, charge and other
parameters. Furthermore, the corresponding thermodynamic entropy is
$S=\mathcal{A}k_{\rm B}c^{3}/4\hbar G$ with $\mathcal{A}$ the horizon area.
The four laws of black hole thermodynamics were subsequently established
Hawking ; Bekenstein ; Bardeen  and it is now extensively believed that black
holes are thermodynamic systems. In contrast to an ordinary fluid system,
black hole entropy is proportional to the horizon area of the black hole
rather than its volume. This brings new peculiar properties to black holes.

An outstanding challenge is that of understanding the microstructure of a
black hole that underlies this thermodynamics. Different theories and models
have been proposed, including string theory Vafa ; Maldacena ; Callan ;
Horowitz , fuzzball models Lunin ; Mathur , and pierced horizons
Rovelli:1996dv , each with its own limitations.

Over the past two decades there has been considerable interest in black hole
phase transitions in anti-de Sitter space Chamblin ; Chamblin2 ; Cognola ;
Kubiznak . Motivated by the interpretation of the cosmological constant as a
thermodynamic pressure Creighton:1995au ; Kastor , four-dimensional charged
AdS black holes were found to exhibit a small-large black hole phase
transition, fully analogous to the liquid-gas phase transition of a Van der
Waals (VdW) fluid Kubiznak . Later studies also showed that this small-large
black hole phase transition is ubiquitous, taking place in most black hole
systems Kubiznak:2016qmn .

It is well known that phase transitions result from the coexistence and
competition among different constituents of a physical system. Phase
transitions thereby provide insight into the properties of a given system’s
microstructure. An effective method, Ruppeiner geometry Ruppeiner , has been
developed to explore this idea. It has been successfully applied to many fluid
systems, including VdW fluids, and many interesting microscopic properties
have been uncovered, shedding new light on thermodynamic systems. Applying
this approach, in recent work Weiw ; Weiwa2 ; WeiWeiWei  we proposed that the
constituent degrees of freedom of black holes have a molecular character.
Specifically, we found that under certain circumstances repulsive
microstructure interactions could exist for black holes, in strong contrast to
VdW fluids Weiwa2 ; WeiWeiWei . However the origin of the repulsive/attractive
interactions remains to be solved.

The present work focuses on this significant problem. We first provide a
microscopic interpretation of Ruppeiner geometry, and then apply it to black
hole systems. We also find the explicit potentials these underlying black hole
molecules experience. We take the following units $G=c=k_{B}=l_{P}=1$.

##  II Microscopic interpretation of Ruppeiner geometry

Ruppeiner observed for various fluid systems that a negative/positive
curvature scalar of the geometry was related to a respective
attractive/repulsive interaction of the fluid molecules Ruppeiner . Employing
this property, Ruppeiner geometry has been used to examine the microstructure
of some complex fluid systems and black holes. Why this empirical observation
works is still a mystery. Since Ruppeiner geometry is closely related to phase
transitions, we begin by providing a microscopic interpretation of Ruppeiner
geometry for the VdW fluid model.

Comparing to an ideal fluid, one needs to take into account the size of and
interaction between the molecules for the VdW fluid, which are described by
the parameters $b$ and $a$, respectively. The equation of state reads

$P=\frac{NT}{V-Nb}-\frac{N^{2}a}{V^{2}}=\frac{T}{v-b}-\frac{a}{v^{2}},$ (2)

where $P$, $T$, $V$, and $v$ are the pressure, temperature, total volume, and
specific volume of the system, and $N$ corresponds to the number of the
molecules contained in the fluid. The equation of state shows a first-order
phase transition, and admits a second order phase transition point at the
critical point ($P_{\rm c}$, $T_{\rm c}$, $v_{\rm c}$)=($\frac{a}{27b^{2}}$,
$\frac{8a}{27b}$, $3b$). Such phenomena can be understood from a microscopic
viewpoint. As is well known, the interaction between two fluid molecules is
given by the Lennard-Jones potential Jones

$\phi=4\phi_{0}\left[\left(\frac{r_{0}}{r}\right)^{12}-\left(\frac{r_{0}}{r}\right)^{6}\right]$
(3)

illustrated in Fig. 1. The red dashed curve and blue solid curve are for the
short-range repulsive interaction and longer-range attractive interaction,
respectively. At $r=r_{0}$, we have $\phi$=0, with the minimum of the
potential $-\phi_{0}$ occurring at $r_{\rm min}=2^{\frac{1}{6}}r_{0}$. A VdW
fluid is described by a “hard-core” model, meaning that the interatomic
distance between molecules cannot be smaller than a molecular diameter
$d=r_{\rm min}$. Hence there is a cutoff at the location of the well, and the
short-range repulsive interaction denoted by the red dashed curve is excluded.
Furthermore, taking the mean field approximation, the characteristic
parameters $a$ and $b$ in the equation of state (2) have a microscopic
interpretation Johnstona

$a=\frac{20\pi r_{0}^{3}\phi_{0}}{9\sqrt{2}},\quad b=\sqrt{2}r_{0}^{3}.$ (4)

where $a>0$, indicating the attractive interaction dominates the fluid
microstructure.

Figure 1: The Lennard-Jones potential. The red dashed curve and blue solid
curve are for the short-range repulsive interaction and longer-range
attractive interaction, respectively.

We now clearly address two issues for a VdW fluid that are often neglected.
The first one is the volume $v\geq b$. When $v=b$, it means the molecules
completely occupy the space of the system. No spatial degrees of freedom are
left, and thus one may conjecture there is no interaction between these
molecules. Furthermore, we can expect that when this condition is violated
repulsive interactions will be present. Second, when $T<\frac{27}{32}T_{\rm
c}=\frac{a}{4b}$, the pressure is negative for a given range of $v$, which may
lead to the repulsive interaction.

Making use of fluctuation theory, Ruppeiner geometry has been used to great
success in revealing the intrinsic microstructure of black holes. In
particular, the normalized Ruppeiner curvature scalar was introduced and shown
for black holes to have universal behaviour near a critical point Weiwa2 ;
WeiWeiWei . Here we consider the relation between the interactions and the
normalized curvature scalar, which for a VdW fluid is given by

$\displaystyle R_{\rm N}=-\frac{2a(v-b)\left(ab+Pv^{3}\right)}{\left(2ab-
av+Pv^{3}\right)^{2}}\;.$ (5)

Obviously, the denominator of $R_{\rm N}$ is non-negative, vanishing when the
system approaches the critical point, which leads to $R_{\rm N}\sim(T_{\rm
c}-T)^{-2}$. Here we are concerned with the numerator of $R_{\rm N}$. Noting
that the parameters $a$ and $b$ are positive, we see that $R_{\rm N}<0$ when
$v\geq b$. This implies an attractive interaction, which is consistent with
the microscopic potential model of the VdW fluid.

We can gain further insight by relaxing the requirements on these parameters.
First, we see that $R_{\rm N}$ vanishes when $a=0$. Since $a$ is a measure of
molecular interaction, vanishing $a$ indicates there is no interaction,
consistent with empirical observations of Ruppeiner geometry. As expected,
from Fig. 1 and Eq. (4), positive/negative $a$ corresponds to a well/barrier,
indicating attractive/repulsive interactions respectively. Since $R_{\rm N}$
has opposite sign to $a$, a negative/positive curvature scalar indicates
respective attractive/repulsive interactions. The absolute value of $R_{\rm
N}$ provides a measure of the strength of these interactions.

Turning to the second term $(v-b)$, if $b=v$, all the space of the fluid
system is occupied by the molecules, so no interactions are allowed. The
corresponding curvature scalar vanishes, consistent with empirical
observations of Ruppeiner geometry. If we do not impose a cutoff at the
location of the well in the Lennard-Jones potential, shifting it to smaller
$r$, then a short-range repulsive interaction is present. There is a range of
negative $(v-b)$, and $R_{\rm N}$ becomes positive as expected.

Finally, considering the term $(ab+Pv^{3})$, it is positive since a fluid
system has positive pressure and volume. But we have noted that when
$T<\frac{a}{4b}$, the pressure becomes negative, yielding positive $R_{\rm
N}$. This case also provides an opportunity for repulsive interactions, and is
a combined effect of both $a$ and $b$.

The VdW parameter space has been shown to have a region where $R_{\rm N}$ is
positive WeiWeiWei , but it completely falls within the coexistence region. It
was therefore excluded since it was not clear that the equation of state was
applicable in the coexistence region. Only the attractive interaction then
dominates in the whole parameter space, which is also consistent with the
microscopic model of a VdW fluid. However if the equation of state is
applicable in this region, the positive region of $R_{\rm N}$ is actually
caused by negative pressure.

In summary, combining with the Lennard-Jones potential, we have given a clear
microscopic interpretation of Ruppeiner geometry for a VdW fluid. We shall now
apply this interpretation to black hole systems, constructing the
corresponding molecular potentials.

##  III Schwarzschild AdS black hole

In an AdS space, Schwarzschild black hole is the most simplest black hole
solution, which is described by the following action

$S=\frac{1}{16\pi G}\int d^{4}x\sqrt{-g}\left(R+\frac{6}{l^{2}}\right),$ (6)

where $l$ is the AdS radius and it is interpreted as the pressure of the black
hole $P=3/(8\pi l^{2})$ Kastor . Solving the corresponding field equations,
the line element reads

$ds^{2}=-f(r)dt^{2}+\frac{1}{f(r)}dr^{2}+r^{2}d\Omega^{2}_{2},$ (7)

where the metric function $f(r)=1-2M/r+r^{2}/l^{2}$ and $M$ denotes the black
hole mass. Requiring the absence of conical singularities at the horizon in
the Euclidean section of the black hole solution, the Hawking temperature can
be easily obtained

$T=\frac{\partial_{r}f(r_{\rm h})}{4\pi}=\frac{1}{4\pi r_{\rm h}}+2Pr_{\rm h}$
(8)

with $r_{\rm h}$ the radius of the black hole horizon, which is actually the
equation of state of the black hole system. It is easy to check that the first
law holds, i.e., $dE=TdS-PdV$, with the respective entropy and thermodynamic
volume $S=\pi r_{\rm h}^{2}$ and $V=\frac{4}{3}\pi r_{\rm h}^{3}$. Noting that
the specific volume $v=2r_{\rm h}l^{2}_{P}$ Kubiznak , we have the number of
the potential black hole molecules

$\displaystyle N=\frac{V}{v}=\frac{2}{3}\pi r_{\rm h}^{2}=\frac{2}{3}S,$ (9)

from which we can write

$\displaystyle P=\frac{NT}{V}-\frac{N^{2}}{2\pi
V^{2}}=\frac{T}{v}-\frac{1}{2\pi v^{2}}$ (10)

for the equation of state for the Schwarzschild AdS black hole. Comparing with
the equation of state of the VdW fluid (2), we can immediately obtain

$\displaystyle a=\frac{1}{2\pi},\quad b=0.$ (11)

After a simple calculation, we give a possible effective potential between two
Schwarzschild AdS black hole molecules

$\displaystyle\phi=\frac{495}{32\sqrt{2}\pi^{2}r_{0}^{3}}\left[\frac{1}{\left(r/r_{0}+\sqrt[6]{2}\right)^{12}}-\frac{1}{\left(r/r_{0}+\sqrt[6]{2}\right)^{6}}\right],$
(12)

which is plotted in Fig. 2. Here the parameter $r_{0}$ has the dimension of
length. Clearly, this gives an attractive interaction for $r/r_{0}\geq 0$. At
$r/r_{0}=0$, one has minimum $\phi_{\rm min}\approx-0.28/r_{0}^{3}$. This
pattern is consistent with that of the normalized scalar curvature given in
Ref. Mannw , where only the attractive interaction dominates for the
Schwarzschild AdS black hole.

Figure 2: The molecular potential (12) for the Schwarzschild AdS black hole.

##  IV Charged AdS black hole

Four-dimensional charged AdS black holes were the first ones found to possess
a small-large black hole phase transition Kubiznak , reminiscent of the
liquid-gas phase transition of a VdW fluid. The line element is still of the
form (7), but with metric function $f=1-2M/r+Q^{2}/r^{2}+r^{2}/l^{2}$, with
$Q$ being the black hole charge. Compared with the Schwarzschild AdS black
hole, the existence of the phase transition is mainly caused by the present of
the black hole charge. The equation of state is

$\displaystyle P=\frac{NT}{V}-\frac{N^{2}}{2\pi V^{2}}+\frac{2N^{4}Q^{2}}{\pi
V^{4}}=\frac{T}{v}-\frac{1}{2\pi v^{2}}+\frac{2Q^{2}}{\pi v^{4}}.$ (13)

A phase transition exists caused by the charge. Regarding the black hole
system as being composed of underlying molecules Weiw ; Weiwa2 , the charge
must be an intrinsic property of these molecules. It is then natural to define
an effective specific charge

$\displaystyle q=N*\frac{Q}{V},$ (14)

from which the equation of state (13) can be cast into a new form

$\displaystyle P=\frac{T}{v}-\frac{1}{2\pi v^{2}}(1-4q^{2}).$ (15)

Comparing with (2) of the VdW fluid, we obtain

$\displaystyle a=\frac{1}{2\pi}(1-4q^{2}),\quad b=0.$ (16)

When the black hole charge vanishes, $q$=0, reducing to (11) for the
Schwarzschild AdS black hole case.

The charged AdS black hole microstructure can be dominated by repulsive
interactions for small black holes Weiwa2 . Following the microscopic
interpretation of Ruppeiner geometry, this can be understood in terms of a
Lennard-Jones like potential that is not cut off at the well. A short-range
repulsive interaction is present that will lead to a dominant repulsive
interaction after a mean-field approximation. When the charge is included, one
can construct the corresponding potential from (12). However as we know, the
molecular potential is not unique for the same equation of state. Here we
would like to model the potential with a new interesting form

$\displaystyle\phi=4\phi_{0}\left[\left(\frac{r_{0}}{r+r_{0}}\right)^{12}-\left(\frac{r_{0}}{r+r_{0}}\right)^{6}\right],$
(17)

where the parameter $r_{0}$ is

$\displaystyle r_{0}^{3}=\frac{495\left(1-4q^{2}\right)}{248\pi^{2}\phi_{0}},$
(18)

which clearly depends on the black hole charge. For this case, the well is
located at $r_{\rm min}=(2^{1/6}-1)r_{0}$. When $0<r<r_{\rm min}$, the slope
of the potential is negative indicating the short-range repulsive interaction.
Otherwise, the interaction is attractive. When $r=r_{\rm min}$, the
interaction vanishes. If we expect it is related with that zero point of the
normalized curvature scalar, we can fix the parameter $\phi_{0}$. From Refs.
Weiwa2 ; WeiWeiWei , we find that when $v=2Q(8\sqrt{2}\pi/3)^{1/3}$, the
curvature scalar vanishes. After setting it to the volume $4\pi r_{\rm
min}^{3}/3$ characterized by the well of the potential, we get

$\displaystyle\phi_{0}\approx\frac{0.0003388}{Q}.$ (19)

We describe the potential $\phi$ (17) for the charged AdS black hole in Fig.
3. Obviously, the short-range repulsive interaction is present at small
$r/Q^{1/3}$. It is worth to point out that the small charge limit is
$q\rightarrow 0$, or equivalently $r_{\rm h}^{2}Q\rightarrow 0$.

Figure 3: The molecular potential (17) for the charged AdS black hole after
fitting the local well with the zero point of the normalized curvature scalar
Weiwa2 ; WeiWeiWei .

##  V Summary

By combining empirical observations Ruppeiner geometry with the microscopic
model of a VdW fluid with the Lennard-Jones potential, we have obtained the
first microscopic interpretation for the constituent degrees of freedom of a
(charged) AdS black hole. We have given for the first time explicit molecular
potentials for these constituents, successfully implementing the potential
black hole molecule model proposed previously Weiw ; Weiwa2 . Our results
address the outstanding problem of why attractive/repulsive interactions
between the microscopic constituents of a thermodynamic system empirically
corresponding negative/positive Ruppeiner curvature scalars.

Our results suggest that the repulsive interaction comes from two aspects. One
is the finite size of the molecular constituents; the other is the cutoff of
the well, where the short-range repulsive interaction should be present in the
potential. In particular, we found that a repulsive interaction is present for
the charged AdS black hole whilst absent for the Schwarzschild AdS black hole.
Note that the molecules have a finite actual size after an approximate
comparison with the VdW fluid Miao , whereas in our case, these black hole
molecules are point particles that may have no actual size but have a specific
volume $v=V/N$, similar to ideal gas molecules.

We close by noting that although our microscopic interpretation of Ruppeiner
geometry was obtained via an AdS black hole, it is applicable to other
asymptotically flat black holes due to the fact that the curvature scalar is
invariant under the coordinate transformation in the thermodynamic parameter
space Ruppeiner . It would be interesting to calculate the molecular
potentials for other black holes with more interesting black hole phase
transition. Moreover, adopting these molecular potentials, it becomes possible
to calculate correlation lengths, which are useful for studying critical
phenomena in black hole phase transitions. These results should provide
considerably more detailed information concerning the underlying degrees of
freedom of black holes.

## Acknowledgements

The authors would like to thank Prof. Hong Zhao for useful discussions. As
this paper was being completed, we became aware of Ref. Dutta , in which the
configuration integral was used to propose an interaction between the black
hole molecules for a dyonic black hole. This work was supported by the
National Natural Science Foundation of China (Grants No. 12075103, No.
11675064, No. 11875151, and No. 12047501), the 111 Project (Grant No. B20063),
and the Natural Sciences and Engineering Research Council of Canada.

## References

  * (1) S. W. Hawking, Particle creation by black holes, Commun. Math. Phys. 43, 199 (1975).
  * (2) J. D. Bekenstein, Black holes and entropy, Phys. Rev. D 7, 2333 (1973).
  * (3) J. M. Bardeen, B. Carter, and S. Hawking, The four laws of black hole mechanics, Commun. Math. Phys.31, 161 (1973).
  * (4) A. Strominger and C. Vafa, Microscopic origin of the Bekenstein-Hawking entropy, Phys. Lett. B 379, 99 (1996), [arXiv:hep-th/9601029].
  * (5) J. M. Maldacena and A. Strominger, Statistical entropy of four-dimensional extremal black holes, Phys. Rev. Lett. 77, 428 (1996), [arXiv:hep-th/9603060].
  * (6) C. G. Callan and J. M. Maldacena, D-brane approach to black hole quantum mechanics, Nucl. Phys. B 472, 591 (1996), [arXiv:hep-th/9602043].
  * (7) G. T. Horowitz and A. Strominger, Counting states of near extremal black holes, Phys. Rev. Lett. 77, 2368 (1996), [arXiv:hep-th/9602051].
  * (8) O. Lunin, S. D. Mathur, AdS/CFT duality and the black hole information paradox, Nucl. Phys. B 623, 342 (2002), [arXiv:hep-th/0109154].
  * (9) O. Lunin, S. D. Mathur, Statistical interpretation of Bekenstein entropy for systems with a stretched horizon, Phys. Rev. Lett. 88, 211303 (2002), [arXiv:hep-th/0202072].
  * (10) C. Rovelli, Black hole entropy from loop quantum gravity, Phys. Rev. Lett. 77, 3288 (1996), [arXiv:gr-qc/9603063 [gr-qc]].
  * (11) A. Chamblin, R. Emparan, C. V. Johnson, and R. C. Myers, Charged AdS black holes and catastrophic holography, Phys. Rev. D 60, 064018 (1999), [arXiv:hep-th/9902170].
  * (12) A. Chamblin, R. Emparan, C. V. Johnson, and R. C. Myers, Holography, thermodynamics and fluctuations of charged AdS black holes, Phys. Rev. D 60, 104026 (1999), [arXiv:hep-th/9904197].
  * (13) M. M. Caldarelli, G. Cognola, and D. Klemm, Thermodynamics of Kerr-Newman-AdS Black Holes and Conformal Field Theories, Class. Quant. Grav. 17, 399 (2000), [arXiv:hep-th/9908022].
  * (14) D. Kubiznak and R. B. Mann, $P$-$V$ criticality of charged AdS black holes, J. High Energy Phys. 1207, 033 (2012), [arXiv:1205.0559[hep-th]].
  * (15) J. D. E. Creighton and R. B. Mann, Quasilocal thermodynamics of dilaton gravity coupled to gauge fields, Phys. Rev. D 52, 4569-4587 (1995), [arXiv:gr-qc/9505007 [gr-qc]].
  * (16) D. Kastor, S. Ray, and J. Traschen, Enthalpy and the mechanics of AdS black holes, Class. Quant. Grav. 26, 195011 (2009), [arXiv:0904.2765[hep-th]].
  * (17) D. Kubiznak, R. B. Mann and M. Teo, Black hole chemistry: thermodynamics with Lambda, Class. Quant. Grav. 34, 063001 (2017), [arXiv:1608.06147 [hep-th]].
  * (18) G. Ruppeiner, Riemannian geometry in thermodynamic fluctuation theory, Rev. Mod. Phys. 67, 605 (1995); erratum ibid 68, 313 (1996).
  * (19) S.-W. Wei and Y.-X. Liu, Insight into the microscopic structure of an AdS black hole from thermodynamical phase transition, Phys. Rev. Lett. 115, 111302 (2015), [arXiv:1502.00386[gr-qc]]; S.-W. Wei and Y.-X. Liu, Erratum: Phys. Rev. Lett. 116, 169903(E) (2016).
  * (20) S.-W. Wei, Y.-X. Liu, and R. B. Mann, Repulsive Interactions and Universal Properties of Charged AdS Black Hole Microstructures, Phys. Rev. Lett. 123, 071103 (2019), [arXiv:1906.10840[gr-qc]].
  * (21) S.-W. Wei, Y.-X. Liu, and R. B. Mann, Ruppeiner Geometry, Phase Transitions, and the Microstructure of Charged AdS Black Holes, Phys. Rev. D 100, 124033 (2019) [arXiv:1909.03887[gr-qc]].
  * (22) S. Dutta, and G. S. Punia, On the Interaction Between AdS Black Hole Molecules, [arXiv:2108.06135 [hep-th]].
  * (23) J. E. Jones, On the determination of molecular fields. -II. From the equation of state of a gas, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character. 106, 463 (1924).
  * (24) D. C. Johnston, Thermodynamic Properties of the van der Waals Fluid, [arXiv:1402.1205 [cond-mat.soft]].
  * (25) S.-W. Wei, Y.-X. Liu, and R. B. Mann, Novel dual relation and constant in Hawking-Page phase transition, Phys. Rev. D 102, 104011 (2020), [arXiv:2006.11503[gr-qc]].
  * (26) Y.-G. Miao and Z.-M. Xu, Thermal molecular potential among micromolecules in charged AdS black holes, Phys. Rev. D 98, 044001 (2018), [arXiv:1712.00545 [hep-th]].