# Understanding Higher-order Structures in Evolving Graphs:
A Simplicial Complex based Kernel Estimation Approach

Manohar Kaul and Masaaki Imaizumi  1Indian Institute of Technology Hyderabad /
2The University of Tokyo

###### Abstract.

Dynamic graphs are rife with higher-order interactions, such as co-authorship
relationships and protein-protein interactions in biological networks, that
naturally arise between more than two nodes at once. In spite of the
ubiquitous presence of such higher-order interactions, limited attention has
been paid to the higher-order counterpart of the popular pairwise link
prediction problem. Existing higher-order structure prediction methods are
mostly based on heuristic feature extraction procedures, which work well in
practice but lack theoretical guarantees. Such heuristics are primarily
focused on predicting links in a static snapshot of the graph. Moreover, these
heuristic-based methods fail to effectively utilize and benefit from the
knowledge of latent substructures already present within the higher-order
structures. In this paper, we overcome these obstacles by capturing higher-
order interactions succinctly as simplices, model their neighborhood by face-
vectors, and develop a nonparametric kernel estimator for simplices that views
the evolving graph from the perspective of a time process (i.e., a sequence of
graph snapshots). Our method substantially outperforms several baseline
higher-order prediction methods. As a theoretical achievement, we prove the
consistency and asymptotic normality in terms of the Wasserstein distance of
our estimator using Stein’s method.

##  1\. Introduction

Numerous types of networks like social (Liben-Nowell and Kleinberg, 2007a),
biological (Airoldi et al., 2006), and chemical reaction networks
(Wegscheider, 1911) are highly dynamic, as they evolve and grow rapidly via
the appearance of new interactions, represented as the introduction of new
links / edges between the nodes of a network. Identifying the underlying
mechanisms by which such networks evolve over time is a fundamental question
that is not yet fully understood. Typically, insight into the temporal
evolution of networks has been obtained via a classical inferential problem
called _link prediction_ , where given a snapshot of the network at time $t$
along with its linkage pattern, the task is to assess whether a pair of nodes
will be linked at a later time $t^{\prime}>t$.

While inferring pairwise links is an important problem, it is oftentimes
observed that most of the real-world graphs exhibit _higher-order group-wise
interactions_ that involve more than two nodes at once. Examples illustrating
human group behavior involve a co-author relationship on a single paper and a
network of e-mails to multiple recipients. In nature too, one can observe
several proteins interacting together in a biological network simultaneously.

In spite of their significance, in comparison to single edge inference,
relatively fewer works have studied the problem of predicting higher-order
group-wise interactions. Benson et al. (2018) originally introduced a
_simplex_ to model group-wise interactions between nodes in a graph. They
proposed predicting a _simplicial closure event_ , whereby an _open simplex_
(with just pairwise interactions between member vertices) transitions to a
_closed simplex_ (where all member vertices participate in the higher-order
relationship simultaneously), in the near future. Figure 2 (Middle) shows an
example of such a transition from an _open triangle_ to a _closed_ one.
Recently, several works have proposed modeling higher-order interactions as
_hyperedges_ in a hypergraph (Xu et al., 2013; Zhang et al., 2018; Yoon et
al., 2020; Patil et al., 2020). Given a hyperedge $h_{t}$ at time $t$, the
inference task is to predict the future arrival of a new hyperedge
$h_{t^{\prime}}$ at time $t^{\prime}>t$, which covers a larger set of vertices
than $h_{t}$ and contains all the vertices in $h_{t}$. Figure 2 (Right)
illustrates this hyperedge prediction task.

Although prediction models based on either _simplicial closure event_
prediction or _hyperedge_ arrival, deal with higher-order structures, they
both fail to capture the _highly complex_ and _non-linear evolution_ of
higher-order structures over time. Both these kinds of models have
limitations. First, they predict structures from a _single static snapshot_ of
the graph, thus not viewing the evolution process of adding new edges as a
time process. Second, their feature extraction is mostly based on popular
heuristics (Adamic and Adar, 2003; Brin and Page, 2012; Jeh and Widom, 2002;
Zhou et al., 2009a; Barabási and Albert, 1999; Bhatia et al., 2019) that work
well in practice but are not accompanied by strong theoretical guarantees. In
addition to the aforementioned shortcomings, hypergraph based methods model
higher-order structures as hyperedges, which omit _lower-dimensional
substructures_ present within a single hyperedge. As a consequence, they
cannot distinguish between various substructure relationships. For example,
hyperedge $[A,B,C]$ in Figure 2 (Right) cannot distinguish between group
relationships like $[[A,B],[B,C],[A,C]]$ (a set of pairwise interactions)
versus $[A,B,C]$ (all $A$, $B$ and $C$ simultaneously in a relationship).
Further, it is important to note that the hypergraph-based approach can be
computationally inefficient. To model a single simplex with $n$ vertices
(inclusive of all possible subsets), the hypergraph model requires $2^{n}$
_explicit_ hyperedges (exponential in $n$), which makes it computationally
prohibitive.

Our Approach: To address the aforementioned problems, we develop a higher-
order structure prediction method, by introducing the following two
techniques: (i) a simplicial complex representation for various events in a
graph, and (ii) a kernel modelling method for prediction with time evolving
events.

A finite collection of non-empty sets $\Delta$ is called a _abstract
simplicial complex_ if, for every set $S\in\Delta$, all its non-empty subsets
$X\subseteq S$ also belong to $\Delta$. The set $S$ is termed a _simplex_ of
$\Delta$. Figure 1 illustrates each simplex comprising of all its lower-
dimensional sub-simplices. We apply this notion to a graph and develop a graph
simplicial complex (GSC). Since the GSC is closed under taking subsets, it
expresses higher-order relationships in a graph in a succinct manner.
Consequently, the _outputs/artifacts_ of a higher-order relationship can be
captured using a _hypergraph_ , but the finer details of _“who interacted with
whom”_ is best captured by a _simplicial complex_.

Figure 1. Illustration of simplices of at most $2$ dimensions. Each simplex
contains all simplices of lower dimensions.

The kernel modeling prediction approach is a method of constructing an
estimator with a kernel function as a time series analysis of events. At the
outset, we regard the _evolving graph_ 111We handle the _incremental model_
(edge insertions only) as opposed to the harder _fully dynamic model_ (edge
insertions and deletions allowed) for which most previous methods too cannot
provide theoretical guarantees. as a _time process_ under the framework of
nonparametric time series prediction. Further, we design a kernel function
based on the features of a GSC to predict the evolution of a given simplex to
a higher-dimensional simplex at a future timestep. To design the kernel
function, we utilize the combination of a _face-vector_ (Björner and Kalai,
2006) (a well-established vector signature in combinatorial topology
literature) and a novel _scoring function_ , which infers the affinity of sub-
simplices based on their past interactions.

Figure 2. [Left] Given a 4-node graph, at time $t$, the 2-simplex $[A,B,C]$
also contains 1-simplices $[A,B],[B,C]$ and $[A,C]$. At time $t^{\prime}>t$,
the 2-simplex _evolves_ (by connecting with $D$) to a 3-simplex $[A,B,C,D]$
which additionally contains 1-simplices $[A,D],[B,D]$ and $[C,D]$, along with
2-simplices $[A,B,D],[A,C,D]$ and $[B,C,D]$. [Middle] Simplex setting with
Benson et al. (2018). The method predicts $[A,B,C]$ (_closed triangle_) at
time $t^{\prime}>t$ from an $[A,B],[B,C]$ and $[A,C]$ (_open triangle_) at
time $t$. [Right] Hypergraph represents $[A,B,C]$ as a hyperedge, without any
of its subsets. It cannot distinguish between $[[A,B],[B,C],[A,C]]$ and
$[A,B,C]$.

Our prediction method has the following advantages.

  * •

It extracts appropriate features of higher-order structures with reasonable
computational complexity.

  * •

Consistency and asymptotic normality of our estimator are proved, while many
of the existing studies do not exhibit such theoretical guarantees.

  * •

In our experiments, our method significantly gains in prediction accuracy and
computational time in comparison to the baselines.

Real-world examples using simplicial complexes: We further motivate our
problem by describing a few significant real-world applications of simplicial
complexes.

(i) Numerical simulation in plasma physics: The _particle-in-cell_ (PIC)
method Evans and Harlow (1957) is a numerical simulation in intense laser-
plasma physics, which traces particles employing $n$-body methods. The
behavior of plasma is largely determined by _particle-particle_ (PP)
interactions Sho (2005). More specifically, group-wise $n$-ary interactions
are studied under several _interaction cross-sections_ Martinez et al. (2019),
where each order-$n$ cross section focuses on a $n$-ary interaction. Since it
is important to distinguish pairwise-particle interactions and a group-wise
interaction between $n$ particles, an approach with simplices is necessary.
(ii) Protein dynamics in structural biology and biochemistry: Proteins are
composed of amino acids linked by covalent peptide bonds, and they are modeled
by _protein structure networks_ Greene (2012), whose structure is directly
correlated to the function of the protein Mannige (2014). Especially, higher-
order interplay between _amino acid groups_ are modeled by “protein sectors”
Halabi et al. (2009). Simplices are suitable for modeling the complicated
interactions in the sectors. A more detailed description of these two examples
is given in the supplementary material (SM).

###  1.1. Related Studies

Single link prediction: Most literature that predicts a single edge/link can
be broadly classified as based on: (i) heuristics, (ii) random-walks, or (iii)
graph neural networks (GNNs). (i) Heuristic methods comprise of Common
neighbors, Adamic-adar (Adamic and Adar, 2003), PageRank (Brin and Page,
2012), SimRank (Jeh and Widom, 2002), resource allocation (Zhou et al.,
2009a), preferential attachment (Barabási and Albert, 1999), persistence
homology based ranking (Bhatia et al., 2019), and similarity-based methods
(Liben-Nowell and Kleinberg, 2007b; Lü and Zhou, 2011). (ii) Random walk based
methods consist of DeepWalk (Perozzi et al., 2014) , Node2Vec (Grover and
Leskovec, 2016a) and SpectralWalk (Sharma et al., 2020). (iii) Finally, for
both link prediction and node classification tasks, recent works are mainly
GNN-based methods such as VGAE (Kipf and Welling, 2016), WYS (Abu-El-Haija et
al., 2018), and SEAL (Zhang and Chen, 2018a).

Higher-order link prediction: Figure 2 summarizes our approach and the related
studies.

Benson et al. (2018) are the first to introduce a _higher-order link
prediction problem_ where they study the likelihoods of future higher-order
group interactions as _simplicial closure_ events (explained earlier). Despite
the novelty, the task proposed by Benson et al. (2018) is more limited than
ours. Our problem setting requires just a single simplex $\sigma$ in order to
predict a higher-dimensional simplex $\tau$, which contains $\sigma$ as a
_face / subset_ , whereas Benson et al. (2018) requires the presence of _all_
constituent $\sigma$ faces in order to predict $\tau$. For example, in Benson
et al. (2018) (also shown in Figure 2 (Middle)) all faces $[A,B]$, $[B,C]$ and
$[A,C]$ (open triangle) need to be present in order to predict a closed
triangle $[A,B,C]$. Contrastingly, in our approach, just a single face/edge
like $[A,B]$ or $[B,C]$ or $[A,C]$, suffices to predict its evolution to
$[A,B,C]$. Figure 2 (Left) illustrates an additional example of our proposal
to predict a $3$-simplex $[A,B,C,D]$ given only one of its faces $[A,B,C]$.

Furthermore, there are studies using _hypergraphs_ which also help naturally
represent group relations (Xu et al., 2013; Zhang et al., 2018; Yoon et al.,
2020; Patil et al., 2020). Especially, to represent higher-order
relationships, Yoon et al. (2020) proposed $n$-projected graphs. For larger
$n$, i.e., higher-order groups, the enumeration of subsets, and keeping track
of node co-occurrences quickly becomes infeasible. In comparison to a
hypergraph, our GSC is closed under taking subsets, which enables us to better
encode more information for improved inference.

##  2\. Preliminary: Graph Simplicial Complex

Figure 3. Example of evolution of GSC $\mathcal{G}$. The yellow triangles are
2-simplices. In the $k$-ball around $1$-simplex $[9,10]$ (in red), a
$1$-simplex [10,7] is added at time $t$.

We start with a general notion of an abstract simplicial complex (ASC), then
define a simplex using ASCs. We specialize this definition to graphs and
define a _graph simplicial complex_ (GSC).

######  Definition 1 (Abstract simplicial complex and simplex).

An abstract simplicial complex (ASC) is a collection $A$ of finite non-empty
sets, such that if $\sigma$ is an element of $A$, then so is every non-empty
subset of $\sigma$. The element $\sigma$ of $A$ is called a simplex of $A$;
its _dimension_ is one less than the number of its elements.

Now, we analyze graphs using the definition of ASCs. Let $G=(V,E)$ be a finite
graph with vertex set $V$ and edge set $E$. A _graph simpicial complex (GSC)_
$\mathcal{G}$ on $G$ is an ASC consisting of subsets of $V$. In particular,
$\mathcal{G}$ is a _collection of subgraphs of G_. With graphs, we denote a
$d$-dimensional simplex (or $d$-simplex) of a GSC by
$\sigma^{(d)}=[v_{0},v_{1},\dots,v_{d}]$. Each non-empty subset of
$\sigma^{(d)}$ is called a _face_ of $\sigma^{(d)}$.

We define several notions related to GSCs, that are useful for describing the
evolution of graphs.

######  Definition 2 (Filtered GSC).

For $I\subset\mathds{N}$, a _filtered GSC_ indexed over $I$ is a family
$(\mathcal{G}_{t})_{t\in I}$ of GSCs such that for every $t\leq t^{\prime}$ in
$I$, $\mathcal{G}_{t}\subset\mathcal{G}_{t^{\prime}}$ holds.

Obviously,
$\mathcal{G}_{t_{0}}\subset\mathcal{G}_{t_{1}}\subset\dots\mathcal{G}_{t_{n}}$
is a discrete filtration induced by the arrival times of the simplices:
$\mathcal{G}_{t_{i}}\setminus\mathcal{G}_{t_{i-1}}=\sigma_{t_{i}}$. This
depicts a higher-order analogue of an evolving graph (incremental model),
which allows _attaching_ new simplices at each time-step to an existing GSC to
build a new GSC. A filtered GSC $\mathcal{G}_{t,p}$ for the last $p$ discrete
time steps is defined as
$\mathcal{G}_{t,p}:=(\mathcal{G}_{t^{\prime}})_{t^{\prime}=t-p}^{t}=(\mathcal{G}_{t-p},\dots,\mathcal{G}_{t})$,
where $\mathcal{G}_{t^{\prime}}\supset\mathcal{G}_{{t^{\prime}}-1}$.

We define a notion for dealing with the neighborhood around a given simplex
$\sigma^{(d)}\in\mathcal{G}$. We introduce a set of all simplices of dimension
$d^{\prime}$ or less from $\mathcal{G}$, i.e.,
$\mathcal{G}^{(d^{\prime})}_{-}:=\\{\sigma^{(d)}\in\mathcal{G}\mid d\leq
d^{\prime}\\}$. $\mathcal{G}^{(0)}_{-}$ is a vertex set and
$\mathcal{G}^{(1)}_{-}$ is the set of edges and vertices. We also write $i\sim
j$ whenever vertices $i$ and $j$ are adjacent in $\mathcal{G}^{(1)}_{-}$, and
write $i\sim_{k}j$ to indicate that vertex $j$ is _$k$ -reachable_ from $i$,
i.e., there exists a path of length at most $k$, connecting $i$ and $j$ in
$\mathcal{G}^{(1)}_{-}$. Then, we define a ball around $\sigma^{(d)}$.

######  Definition 3 ($k$-ball centered at vertex and simplex).

At time $t$, we define a $k$-ball centered at vertex $i$ by
$B_{t,k}(i):=\\{j:i\sim_{k}j\text{ and }i,j\in\mathcal{G}^{(0)}_{t-}\\}$, and
a $k$-ball centered at a simplex $\sigma^{(d)}$ as
$B_{t,k}(\sigma^{(d)}):=\bigcup_{i:\mathrm{Vert}(\sigma^{(d)})}B_{t,k}(i)$,
where $\mathrm{Vert}(\sigma^{(d)})$ denotes all vertexes in $\sigma^{(d)}$.

Now, we define a _sub-complex_
$\mathcal{G}^{\prime}_{t}(\sigma^{(d)})\subseteq\mathcal{G}_{t}$ as the GSC
that contains all the simplices in $\mathcal{G}_{t}$ spanned by the vertices
in the $k$-ball $B_{t,k}(\sigma^{(d)})$.

##  3\. Predicting Higher-order Simplices

We consider the prediction of a simplex’s arrival in the setting described
below. Consider a filtered GSC $\mathcal{G}_{t,p}$. At time $t$, given a
$d$-dimensional simplex $\sigma^{(d)}=[v_{0},\cdots,v_{d}]$, we predict the
formation of a $(d+1)$-simplex
$\tau^{(d+1)}=[v_{0},\cdots,v_{d},\widetilde{v}]$ with a new vertex
$\widetilde{v}\in B_{t,k}(\sigma^{(d)})$. Here, we restrict $\widetilde{v}$ to
be $k$-reachable from $\sigma^{(d)}$. To find out which simplices and vertices
are most likely to appear in $\tau^{(d+1)}$ at time $t+1$, we need to design
features for $\sigma^{(d)}$ and $\widetilde{v}$.

###  3.1. Feature Design for Simplex

We develop a feature design of a simplex associated with the notion of
$k$-balls. The design is organized into two main elements: (i) a face-vector
with a sub-complex, and (ii) a scoring function.

(i) Face vector of sub-complex: We first define a face-vector of a fixed GSC.
The face-vector is an important _topological invariant_ 222A topological
invariant is a property that is preserved by _homeomorphisms_. of the GSC.

######  Definition 4 (face-vector).

A combinatorial statistic of $\mathcal{G}$ is the face-vector (or $f$-vector)
as $f(\mathcal{G})=(f_{-1},f_{0},\dots,f_{d-1})$, where
$f_{k}=f_{k}(\mathcal{G})$ records the number of $k$-dimensional faces
$\sigma^{(k)}\in\mathcal{G}$.

We then define the feature of a simplex $\sigma^{(d)}$ at time $t$, denoted by
$N_{t}(\sigma^{(d)})=f(\mathcal{G}^{\prime}_{t}(\sigma^{(d)}))$. In words, the
feature is compactly represented as the face-vector of sub-complex
$\mathcal{G}^{\prime}_{t}(\sigma^{(d)})$. Our face-vector representation of a
node’s neighborhood can be considered as a _higher-order analogue_ of the
_Weisfeiler-Lehman_ (WL) kernel Shervashidze et al. (2011) on unlabeled
graphs, which for each vertex, iteratively aggregates the vertex degrees of
its immediate neighbors to compute a unique vector of the target vertex that
captures the structure of its extended neighborhood.

(ii) Scoring function: The purpose of this function is to extract the features
of $\widetilde{v}$ using its proximity to $\sigma^{(d)}$. To this end, we
begin by describing _affinity_ between two vertices. Given two vertices
$v,v^{\prime}\in\mathcal{G}^{(0)}_{-}$, we denote by $s(v,v^{\prime})$ the
weighted sum of all past co-occurrences of vertices $v$ and $v^{\prime}$ in
$\sigma^{(d)}$, where the weight is $d$ from $\sigma^{(d)}$. We then devise a
scoring function $h(\cdot,\cdot)$ that assigns an integral score to the
possible introduction of a vertex $\widetilde{v}$ to a $d$-simplex
$\sigma^{(d)}=[v_{0},\cdots,v_{d}]$ as

$\displaystyle
h_{t}(\sigma^{(d)},\widetilde{v})=\sum_{i=0}^{d}s(v_{i},\widetilde{v}).$ (1)

It describes _higher co-occurrence_ between $\sigma^{(d)}$ and $\widetilde{v}$
at time $t$, indicates a higher likelihood of forming a $(d+1)$-simplex
$\tau^{(d+1)}=[\sigma^{(d)},\widetilde{v}]$ together at a future time $t+1$.
We give a higher score to past co-occurrences of vertex pairs in higher
dimensional simplices.

Feature vector: Finally, for a given $d$-simplex $\sigma^{(d)}$ at time $t$,
and a possible introduction of a new vertex $\widetilde{v}\in
B_{t,k}(\sigma^{(d)})$, we assign a feature vector

$\displaystyle
F_{t}(\sigma^{(d)},\widetilde{v})=(N_{t}(\sigma^{(d)}),h_{t}(\sigma^{(d)},\widetilde{v})).$
(2)

We denote the set of all such possible $(\sigma^{(d)},\widetilde{v})$ pairs
with their corresponding feature vectors equal to $F$ as
$P_{t}(\sigma^{(d)},F)$. Furthermore, among the pairs in
$P_{t}(\sigma^{(d)},F)$, we denote by $P^{\tau}_{t}(\sigma^{(d)},F)$ those set
of pairs with feature vectors equal to $F$ that _actually_ form
$\tau^{(d+1)}=[\sigma^{(d)},\widetilde{v}]$ at time $t$, i.e., $\sigma^{(d)}$
appears as a face in $\tau^{(d+1)}$ at time $t$. Note the distinction that not
all $d$-simplices counted in $P_{t}(\sigma^{(d)},F)$ end up being _promoted_
to higher $(d+1)$-simplices in the next time step. Table 1 provides a list of
notations.

Table 1. Notation table  Basic |
---|---
$G=(V,E)$ |  graph with vertex set $V$ and edge set $E$
$\mathcal{G}$ |  graph simplicial complex (collection of subgraphs of $G$)
$\sigma^{(d)}=[v_{0},...,v_{d}]$ |  $d$-dimensional simplex ($d$-simplex)
$\tau^{(d+1)}=[v_{0},...,v_{d},\widetilde{v}]$ |  $d+1$-simplex for prediction
$\mathcal{G}_{t,p}=\\{\mathcal{G}_{t-p},...,\mathcal{G}_{t}\\}$ |  GSCs from the previous $p$ time steps
$\mathcal{G}_{-}^{(d^{\prime})}=\\{\sigma^{(d)}\in\mathcal{G}\mid d\leq d^{\prime}\\}$ |  set of simplices of dimension $d^{\prime}$ or less
Local simplex |
$B_{t,k}(i)=\\{j:i\sim_{k}j,~{}i,j\in\mathcal{G}_{t-}^{(0)}\\}$ |  $k$-ball centered at vertex $i$
$B_{t,k}(\sigma^{(d)})=\bigcup_{i:\mathrm{Vert}(\sigma^{(d)})}B_{t,k}(i)$ |  $k$-ball centered at simplex $\sigma^{(d)}$
$\mathcal{G}_{t}^{\prime}(\sigma^{(d)})$ |  all simplices from $\mathcal{G}_{t}$ spanned by vertices in $B_{t,k}(\sigma^{(d)})$
Feature of simplex |
$f_{k}=f_{k}(\mathcal{G})$ |  total number of $k$-simplices $\sigma^{(k)}\in\mathcal{G}$
$f(\mathcal{G})=(f_{-1},f_{0},...,f_{d-1})$ |  face-vector of $\mathcal{G}$
$N_{t}(\sigma^{(d)})=f(\mathcal{G}_{t}^{\prime}(\sigma^{(d)}))$ |  feature of $\sigma^{(d)}$
$s(v,v^{\prime})$ |  weighted sum of past co-occurrences of $v,v^{\prime}\in\sigma^{(d)}$
$h_{t}(\sigma^{(d)},v)=\sum_{i=0}^{d}s(v_{i},\widetilde{v})$ | scoring function
$F_{t}(\sigma^{(d)},\widetilde{v})=(N_{t}(\sigma^{(d)}),h_{t}(\sigma^{(d)},\widetilde{v}))$ | feature vector
$P_{t}(\sigma^{(d)},F)$ |  set of $\sigma^{(d)},\widetilde{v}$ with corresponding feature $F$

###  3.2. Prediction Model and Kernel Estimator

For the prediction, we define an indicator variable that displays the
appearance of a new simplex. Given a $d$-simplex
$\sigma^{(d)}=[v_{0},\cdots,v_{d}]\in\mathcal{G}_{t}$, the arrival at time
$t+1$ of a $(d+1)$-simplex $\tau^{(d+1)}=[v_{0},\cdots,v_{d},\widetilde{v}]$
with a new vertex $\widetilde{v}\in B_{t,k}(\sigma^{(d)})$ is captured by the
following variable

$Y_{t+1}(\tau^{(d+1)}):=\begin{cases}1&\text{if $\sigma^{(d)}$ is a face of
$\tau^{(d+1)}$ },\\\ 0&\text{otherwise}.\end{cases}$ (3)

Prediction model: Our approach for the prediction is to model the indicator
variable. Namely, we assume that the indicator variable follows the following
distribution:

$\displaystyle
Y_{t+1}(\tau^{(d+1)})\mid\mathcal{G}_{t,p}\sim\text{Bernoulli}(g(F_{t}(\sigma^{(d)},\widetilde{v})))$
(4)

where $0\leq g(\cdot)\leq 1$ is a function of the feature vector
$F_{t}(\sigma^{(d)},\widetilde{v})$. In words, the indicator variable
$Y_{t+1}(\tau^{(d+1)})$ is Bernoulli distributed with success probability
given by function $g(F_{t}(\sigma^{(d)},\widetilde{v}))$, conditioned on
having seen the last $p$ states of GSC $\mathcal{G}_{t}$. This model describes
that the appearance probabilities for two simplices $\sigma_{i}$ and
$\sigma_{j}$ are likely to be similar, if their feature vector is also
similar.

Estimator with Kernels: We utilize a kernel method to estimate the success
probability of the model (Equation 4) based on observed simplices at time $t$.
Let $\mathcal{G}_{t}(d)$ be a set of $d$-dimensional simplices at time $t$
from $\mathcal{G}_{t,p}$. Also, for brevity, let $F$ represent a feature
$F_{t}(\sigma^{(d)},\widetilde{v}_{m})$ with some $t,\sigma^{(d)}$ and
$\widetilde{v}_{m}$ subject to $\widetilde{v}_{m}\in
B_{t,k}(\sigma_{i}^{(d)})$. Let $\lVert F-F^{\prime}\rVert_{1}$ denote the
$L_{1}$-distance between two feature vectors, and also define a $L_{1}$-ball
$\Gamma(F,\delta):=\\{F^{\prime}:\|F-F^{\prime}\|_{1}\leq\delta\\}$.

We define our kernel function $K(\cdot,\cdot)$ as follows. With two features
$F$ and $F^{\prime}$, we define it as

$\displaystyle
K\left(F,F^{\prime}\right):=\frac{\mathds{I}\left\\{F=F^{\prime}\right\\}+\beta\mathds{I}\left\\{\lVert
F-F^{\prime}\rVert_{1}\leq\delta\right\\}}{1+\beta|\Gamma(F,\delta)|},$ (5)

where $\beta>0$ is the bandwidth parameter, and $\mathds{I}$ is the indicator
function. Given the feature $F$ with only integer components, we are
interested in only those _close by feature vectors_ that are either exactly
the same as $f$ or lie within an $L_{1}$-ball of radius $\delta$ centered on
$f$. This explains the choice of our kernel function with discrete indicator
variables.

Now, we define our estimator. At time $T$, we fix $\sigma^{(d)}$ and
$\widetilde{v}$ and set $F=F_{T}(\sigma^{(d)},\widetilde{v})$. Then, we
consider a set of observed feature at time $T$ as

$\displaystyle\mathcal{F}_{T}:=\bigcup_{t^{\prime}=T-p}^{T}\left\\{[\sigma^{(d)}_{j},\widetilde{v}_{n}]:\sigma^{(d)}_{j}\in\mathcal{G}_{t^{\prime}}{(d)},\widetilde{v}_{n}\in
B_{t^{\prime},k}(\sigma_{j}^{(d)})\right\\}.$

Our estimator of $g(F)$ in (Equation 4) is written as follows:

$\displaystyle\widetilde{g}_{T}\left(F\right):=\frac{\sum_{F^{\prime}\in\mathcal{F}_{T}}K\left(F,F^{\prime}\right)\cdot
Y_{t^{\prime}+1}(F^{\prime})}{\sum_{F^{\prime}\in\mathcal{F}_{T}}K\left(F,F^{\prime}\right)}.$
(6)

Time-complexity of our estimator: As both the face-vectors and counts
$P_{t+1}$ / $P^{(\tau)}_{t+1}$ (updated in a data cube) are computed
simultaneously, they incur the same time overhead. At time $t$, for a
$d$-simplex $\sigma^{(d)}$ it takes $O(|V_{k}|+|E_{k}|)$ time to compute a
$k$-ball around $\sigma^{(d)}$, where $V_{k}$ and $E_{k}$ are the set of
vertices and edges in the $k$-hop subgraph of a vertex. We must check this
$k$-ball of $\sigma^{(d)}$ against $p|\mathcal{G}_{t-}^{(d)}|$ number of
simplices (with dimension at most $d$) from the previous $p$ time steps, by
intersecting against them to get counts for the face vector and data cube.
Each intersection test takes $O(|V_{k}|)$ time. Recall that
$f_{d}(\mathcal{G}_{t})$ denoted the total number of $d$-simplices in
$\mathcal{G}_{t}$. $V_{k}$ and $E_{k}$ are the set of vertices and edges in
the $k$-hop subgraph of a vertex. We must check this $k$-ball of
$\sigma^{(d)}$ against $p|\mathcal{G}_{t-}^{(d)}|$ number of simplices (with
dimension at most $d$) from the previous $p$ time steps. So, the entire
computation across $T$ time chunks has a time complexity of
$O(Tpf_{d}(\mathcal{G}_{t})(|V_{k}|+|E_{k}|))|\mathcal{G}_{t-}^{(d)}|)$.

Storage-complexity of our estimator: The storage cost of evolving graphs per
window of size $p$ is $O(|\mathcal{G}_{t,p}|)$. For a single simplex, our
estimator requires storing: (i) a pair of integer counts, namely
$(P_{t}(\cdot,\cdot),P^{\tau}_{t}(\cdot,\cdot))$, in a datacube, which costs
$O(1)$ and (ii) a $d+1$-dimensional face vector which takes storage $O(d+1)$.
Recall, the total number of simplices is denoted by $f_{d}(\mathcal{G}_{t})$.
We arrive at a total storage cost of
$O(df_{d}(\mathcal{G}_{t})+|\mathcal{G}_{t,p}|)$.

##  4\. Theoretical Property of the Estimator

We show that our estimator has theoretical validity: (i) consistency and (ii)
asymptotic normality. (i) The consistency guarantees $\widetilde{g}_{T}$
achieves zero error as $T$ increases by converging to $g$. (ii) The asymptotic
normality implies that the error $\widetilde{g}_{T}-g$ converges to a normal
distribution, which is useful to evaluate the size of the error and can be
applied to statistical tests and confidence analysis. Both properties are very
important in statistics (Van der Vaart, 2000).

###  4.1. Consistency

We study the consistency of our estimator. To discuss the property of the
estimators with GSCs, it is necessary to organize the Markov property, that
the GSC evolution process clearly exhibits. It is well-known that there exists
a set of irreducible closed communication classes $C$ in the state space
$\mathcal{S}$. We denote the time of entering class $C$ by $T_{C}$ and the
event as $\mathcal{E}(T_{C})$. Let $S_{C}$ denote the event $S_{t}\in C$,
where $S_{t}$ is the state of the Markov chain at time $t$. Then,
$\mathcal{E}(T_{C})\cap S_{C}$ is the event that the chain enters class $C$ at
time $T_{C}$ and remains in that communication class indefinitely.

With the event $S_{C}$, we provide the bias-variance decomposition of
$\widetilde{g}_{T}-g$, which is common for theoretical analysis of estimators.
The bias represents an error due to the expressive power of the model and the
variance represents the over-fitting error due to algorithm uncertainty. By
analyzing these terms separately, we can analyze the overall prediction error.
We define two functions as

$\displaystyle\widehat{h}_{T}(F)=(T-p)^{-1}\sum_{t=p}^{T-1}\sum_{j=1}^{\lvert\mathcal{G}_{t}(d)\rvert}{\lvert\mathcal{G}_{t}(d)\rvert}^{-1}{\lvert
P^{\tau}_{t+1}(\sigma^{(d)}_{j},F)\rvert},$
$\displaystyle\widehat{d}_{T}(F)=(T-p)^{-1}\sum_{t=p}^{T-1}\sum_{j=1}^{\lvert\mathcal{G}_{t}(d)\rvert}{\lvert\mathcal{G}_{t}(d)\rvert}^{-1}{\lvert
P_{t+1}(\sigma^{(d)}_{j},F)\rvert}.$

For the sake of brevity, we fix $F$ and denote $\widetilde{g}_{T}(F)$ by
$\widetilde{g}_{T}$. Similarly, $g,\widehat{h}_{T}$ and $\widehat{d}_{T}$ are
used. Also, we define a term $B_{T}(F,C)=\mathds{E}[\widehat{h}_{T}\mid
S_{C}]/\mathds{E}[\widehat{d}_{T}\mid S_{C}]-g$. Then, by Proposition 2 in SM,
we decompose $\widetilde{g}_{T}-g$ as

$\displaystyle\widetilde{g}_{T}-g=\mathcal{V}_{T}+\mathcal{B}_{T},$ (7)

where $\mathcal{V}_{T}$ is a variance term

$\displaystyle\mathcal{V}_{T}:=\\{[\widehat{h}_{T}-g\widehat{d}_{T}]-\mathds{E}[\widehat{h}_{T}-g\widehat{d}_{T}\mid
S_{C}]\\}/{\widehat{d}_{T}},$

and $\mathcal{B}_{T}$ is a bias term

$\displaystyle\mathcal{B}_{T}:=B_{T}(F,C)\mathds{E}[\widehat{d}_{T}\mid
S_{C}]/\widehat{d}_{T}.$

To bound the variance term $\mathcal{V}_{T}$, we make an assumption that our
Markov chain $X_{t}$ exhibits a _$\alpha$ -mixing_ property which describes a
dependent property of the dynamic process. It is one of the most common and
well-used assumptions describing time-dependent processes including dynamic
graphs (Sarkar et al., 2014b). Precisely, we present the definition of
_$\alpha$ -mixing_:

######  Definition 5 ($\alpha$-mixing).

A stochastic process $X_{t}$ is $\alpha$-mixing, if a coefficient $\alpha(r)$,
defined as

$\displaystyle\alpha(r)=\sup_{|t_{1}-t_{2}|\geq r}$
$\displaystyle\\{\lvert\mbox{Pr}(A\cap
B)-\mbox{Pr}(A)\mbox{Pr}(B)\rvert:A\in\Sigma(X^{-}_{t_{1}}),B\in\Sigma(X^{+}_{t_{2}})\\},$

satisfies $\alpha(r)\rightarrow 0$ as $r\to\infty$. Here,
$\Sigma(X^{-}_{t_{1}})$ and $\Sigma(X^{+}_{t_{2}})$ are the sigma-algebras of
past and future events of the stochastic process up to and including $t_{1}$.

This definition implies that time-dependent processes get _close to
independent_ as time passes. That is, events at time $t$ and $t+10000$ are
close to independent, while events at time $t$ and $t+1$ can be correlated.
The simplest example is the evolution of stock prices in financial markets:
the movement of stock prices today does not correlate with the movement of
stock prices $10$ years ago.

To bound the bias term $\mathcal{B}_{T}$, we impose a smoothness condition on
$g$. A similar assumption is often used in the problem of predicting links
(e.g. Assumption 1 in Sarkar et al. (2014a)). Our assumption is a general and
weaker version of the common assumption.

######  Assumption 1 (Smoothness on $g$).

There exists a function $\kappa:\mathbb{R}\to\mathbb{R}$ in the Schwartz space
(i.e. it is infinitely differentiable and converging to zero faster than any
polynomial as $x\to\pm\infty$) with $b>0$ such that
$|g(F)-g(F^{\prime})|=O(\kappa(-\|F-F^{\prime}\|_{1}/b)),~{}\mbox{as}~{}b\to
0,~{}\forall F,F^{\prime}.$

Then, we prove the consistency of $\widetilde{g}_{T}$.

######  Theorem 1 (Consistency).

Suppose that the GSC filtration process is $\alpha$-mixing, $\beta=o(1)$, and
Assumption 1 holds. Then, for any $F$ and conditional on $S_{C}$, our
estimator $\widetilde{g}_{T}(F)$ is well-defined with probability tending to
$1$, and $|\widetilde{g}_{T}(F)-g(F)|\overset{p}{\to}0$ holds as
$T\rightarrow\infty$.

###  4.2. Asymptotic Normality

We show the asymptotic normality of the proposed estimator. That is, we prove
that the error of the estimator converges weakly to a normal distribution.
This property allows for more detailed investigations, such as correcting for
errors in estimators or performing statistical tests.

Technically speaking, we develop a distribution approximation result with
Wasserstein distance (Villani, 2008) and Stein’s method (Stein et al., 1986)
to handle the dependency property of GSCs. We are interested in approximating
a random variable $Z_{n}$ by a Gaussian random variable, where $Z_{n}$ is a
sum of $n$ mean-centered random variables $\\{A_{i}\\}_{i=0}^{n}$, where
$A_{i}$ corresponds to a random variable which depends on the $i$-th
$d$-simplex $\sigma_{i}^{(d)}$ in our GSC, which is _dependent_ on other
$d$-simplices whose neighborhoods largely overlap with that of
$\sigma_{i}^{(d)}$. Here, let $d_{w}$ be the Wasserstein distance between the
underlying distributions of the random variables, and $N$ is a standard
Gaussian variable. Then, we develop the general results for $Z_{n}$. We
provide the following theoretical result. Its formal statement is Theorem 3,
which is deferred to SM due to its complexity.

######  Proposition 1 (Gaussian approximation for dependent variables; Simple
version of Theorem 3).

Suppose $Z_{n}=\sum_{i=0}^{n}A_{i}$, where $\\{A_{i}\\}_{i=0}^{n}$ is
generated by zero-mean random variables $X_{0},X_{1},\cdots,X_{n}$ satisfying
the $\alpha$-mixing condition, such as $A_{i}=X_{i}/B_{n}$ with
$B_{n}=\mathds{E}[\sum_{i=0}^{n}X_{i}^{2}]$. Also, suppose that
$\mbox{Pr}\\{|X_{i}|\leq L\\}=1$ holds for $i=1,\cdots,n$ with some constant
$L>0$. Then, with an existing finite constant $C>0$, we have
$d_{w}(Z_{n},N)\leq
C\left(\sum_{i=1}^{n}\mathds{E}|A_{i}|^{3}+\frac{TL^{3}}{B_{n}^{3}}\sum_{i=1}^{n-1}n\alpha(n)\right)$.

This result extends Sunklodas (2007) in the sense of Markov chains on GSCs
that satisfy the $\alpha$-mixing condition. This extension makes it possible
to study the _contributing effect_ of neighboring $d$-simplices (represented
as a set of weakly dependent r.v.’s) on a central $d$-simplex.

Table 2. AUC scores and runtimes for baselines versus our method’s estimator for $d=1$ and $d=2$. | Enron | Contact | NDC | EU |  | Enron | Contact | NDC | EU
---|---|---|---|---|---|---|---|---|---
| ($d=1$) AUC / runtime (sec) |  | ($d=2$) AUC / runtime (sec)
AA | 0.54 / 0.18 | 0.57 / 0.18 | 0.30 / 0.38 | 0.61 / 0.15 |  | 0.31 / 0.16 | 0.33 / 0.20 | 0.47 / 0.25 | 0.25 / 0.28
JC | 0.42 / 0.20 | 0.63 / 0.2 | 0.16 / 0.37 | 0.65 / 0.20 |  | 0.41 / 0.16 | 0.44 / 0.21 | 0.23 / 0.24 | 0.32 / 0.27
PA | 0.55 / 0.15 | 0.60 / 0.28 | 0.55 / 0.11 | 0.38 / 0.09 |  | 0.52 / 0.15 | 0.63 / 0.20 |  0.74 / 0.24 | 0.34 / 0.18
NV | 0.45 / 74 | 0.25 / 155 | 0.30 / 406 |  0.67 / 280 |  | 0.49 / 71 | 0.45 / 155 | 0.49 / 3989 | 0.40 / 374
SL | 0.54 / 152 |  0.91 / 241 | 0.33 / 260 | 0.57 / 431 |  | 0.48 / 152 | 0.54 / 241 | 0.40 / 260 | 0.29 / 431
TT | 0.62 / 180 | 0.80 / 202 | 0.42 / 285 | 0.61 / 420 |  |  0.70 / 180 |  0.78 / 205 | 0.45 / 285 |  0.67 / 422
TN |  0.67 / 212 | 0.84 / 241 | 0.48 / 280 | 0.55 / 512 |  | 0.62 / 212 | 0.70 / 241 | 0.46 / 281 | 0.64 / 512
HP | 0.26 / 22 | 0.64 / 142 |  0.57 / 54 | 0.45 / 17 |  | 0.26 / 22 | 0.76 / 144 | 0.45 / 58 | 0.41 / 18
Ours |  0.88 / 1.08 |  0.87 / 7.56 |  0.78 / 5.84 |  0.83 / 1.48 |  |  0.94 / 2.25 |  0.83 / 2.76 |  0.96 / 0.76 |  0.80 / 0.642
|  ($\beta$=1) |  ($\beta$=0.1) |  ($\beta$=0.1) |  ($\beta$=10) |  |  ($\beta$=0.01) |  ($\beta$=0.01) |  ($\beta$=0.01) |  ($\beta$=10)

We provide the asymptotic normality of our estimator. It shows that our
estimator converges to a normal distribution in terms of the Wasserstein
distance, which leads to weak convergence. To achieve the result, we utilize
the decomposed terms $\mathcal{V}_{T}$ and $\mathcal{B}_{T}$ from (7). Then,
we regard $\mathcal{V}_{T}$ as a sum of dependent random variables and apply
the result developed in Proposition 1. Let $\sigma_{c}^{2}$ be a limit of
variance of the numerator in $T^{-1/2}\mathcal{V}_{T}$ as $T\to\infty$. We
recall that $S_{C}$ denotes the event $S_{t}\in C$, where $S_{t}$ is the state
of the Markov chain at time $t$.

######  Theorem 2 (Asymptotic Normality).

Suppose that Assumption 1 holds, the GSC filtration process is
$\alpha$-mixing, and $\sigma_{c}>0$. If $\beta=o(T^{-1/2})$ and
$b=o(T^{-1/2})$, then, for any $F$ and conditioned on $S_{C}$, the following
holds:
$\sqrt{T}(\widetilde{g}_{T}(F)-g(F))\overset{d}{\to}\mathcal{N}(0,\sigma_{c}^{2}/R(C)^{2}),\mbox{~{}as~{}}T\to\infty.$

By using this property, we can make detailed inferences based on the
distribution of the estimation error. For example, it is possible to create
confidence intervals for predictions and perform statistical tests to
rigorously test hypotheses about simplex arrivals.

##  5\. Real-World Data Experiments

We empirically evaluate the performance of our proposed estimator on real-
world dynamic graphs compared to baselines. The basic premise in our
experiments is to capture local and higher-order properties surrounding a
$d$-simplex up to time $t$ to predict the appearance of a new $(d+1)$-simplex
at time $t^{\prime}>t$, which contains $\sigma^{(d)}$ as its _face_. Note that
we compare our method to other _closely related_ methods that were designed to
solve different structure prediction tasks.

Datasets: We report results on real-world dynamic graph datasets sourced from
Benson et al. (2018). Each dataset contains $n$ nodes, $m$ formed edges, and
$x$ _timestamped simplices_ (represented as a set of nodes). There are four
datasets named Enron ($n=143$, $m=1.8$K, $x=5$K), EU ($n=998$, $m=29.3$K,
$x=8$K), Contact ($n=327$, $m=5.8$K, $x=10$K), and NDC (National Drug Code)
($n=1.1$K, $m=6.2$K, $x=12$K).

Experimental setup: We first ordered by arrival times and grouped the
timestamped simplices into $T$ time slices. For most of our experiments, $T$
was set to $20$, except for $d=2$, where $T$ was set to $6$ and $12$ for _EU_
and _NDC_ , respectively. Then, we randomly sampled a set of $d$-simplices
from the time slices in the range $[1,T-1]$. Those $d$-simplices paired with a
vertex that successfully formed a face in a $(d+1)$-simplex in the $T$-th time
slice were classified as _positive samples_ , while the rest were deemed as
_negative samples_. We picked an equal number of positive and negative samples
for evaluation. For $K$-fold cross-validation for $\beta$, we swapped the
$T$-th time slice with one of the $K$ slices preceding the $T$-th time slice
for each fold. $K$ was set to $3$. All experiments where repeated $10$ times
and average AUC scores and runtimes are reported.

Compared methods: As naive baselines, we averaged the results of single-edge
prediction methods, where a new edge would form between each node in the
$d$-simplex and the vertex to be paired with. Specifically, we compare our
estimator with: (i) _heuristic_ (_Adamic-Adar_ (AA) (Adamic and Adar, 2001),
_Jaccard Coefficient_ (JC) (Salton and McGill, 1986), and _Preferential
attachment_ (PA) (Mitzenmacher, 2004), (ii) _deep-learning_ based (_Node2vec_
(NV) (Grover and Leskovec, 2016b), and _SEAL_ (SL) (Zhang and Chen, 2018b)),
and (iii) _temporal graph network_ based (_TGAT_ (TT) da Xu et al. (2020) and
_TGN_ (TN) Rossi et al. (2020)) link prediction methods. We note that (Benson
et al., 2018) for predicting a “simplicial closure” has the closest motivation
to our method, yet has divergent objectives, therefore we omit comparison to
their work. For _hyper-edge prediction_ (HP), we picked the recent most
representative work by Yoon et al. (2020) to compare against, although this
work only works for static non-evolving hypergraphs.

###  5.1. Results and Discussion

We averaged the classification accuracy and runtimes of our estimator and the
baselines. We performed two sets of experiments on the arrival of a
$(d+1)$-simplex and summarize it in Table 2 for $d=\\{1,2\\}$. We also report
the bandwidth $\beta$ for our estimator selected by cross-validation.

Predicting $2$-simplex ($d=1$): We observe that our method is nearly two
orders of magnitude faster than the deep learning based methods (NV and SL)
and nearly an order of magnitude faster than the hypergraph prediction method
(HP). While the single edge heuristic methods are relatively faster, their AUC
scores are not comparable to our method’s AUC scores. Also, we achieve nearly
$30\%$ improvement (in Enron) over the next best performing prediction method.

$1$$2$$3$$4$$5$$6$$7$$8$$0.6$$0.7$$0.8$$0.9$$1$Predicted Simplex
DimensionAUCEnronContactNDCEU Figure 4. AUC score predicted by our estimator
for future formation of a $d$-simplex, given a $d-1$-simplex.

Predicting $3$-simplex ($d=2$): The gap in AUC scores between our method and
the baselines are far more pronounced. Our runtimes also improve due to the
far fewer number of simplices with dimensions exceeding $3$. As observed in
Yoon et al. (2020) about slight drops in accuracy for higher-dimensional
hyper-edges, we also note that in HP, the AUC score remains the same or drops
slightly compared to prediction at $d=1$.

Advantage of higher dimensional simplices: We perform additional experiments
by increasing $d$ from $1$ to $8$ and show that handling high-dimensional
simplices exhibits high prediction accuracy. In Figure 4, the prediction is
basically improves as $d$ increases.

Empirical summary: Traditional estimators fail to accurately capture the rich
latent information present in higher-order structures (and their sub-
structures) that evolve over time. Our estimator succinctly captures this
information via the $f$-vector and weighted scoring of $(\sigma^{(d)},v)$ pair
formation depending on the dimension of the simplex in which the pair co-occur
in the past.

##  6\. Conclusion

We modeled the higher-order interaction as a simplex and demonstrated a novel
kernel estimator to solve the higher-order structure prediction problem. From
a theoretical standpoint, we proved the consistency and asymptotic normality
of our estimator. We empirically argue that our estimator outperforms
hypergraph based and higher-order link prediction baselines from both
heuristic and deep-learning based pairwise link prediction methods.

##  Appendix A Further Description on Example of Applications

###  A.1. Plasma Physics

The _particle-in-cell_ (PIC) method Evans and Harlow (1957) is a numerical
simulation to model and study the evolution of the kinetic and collective
phenomena at play in intense laser-plasma physics. It traces the trajectory of
particles employing $n$-body methods, while solving Maxwell’s equations on a
Eulerian grid. The behavior of plasma is largely determined by the interaction
(i.e., both _dynamic collisions_ and _static contact_) between its constituent
particles, generally termed as _particle-particle_ (PP) interactions Sho
(2005). PP interactions are best modeled by weighted dynamic graphs, whose
edge weights represent interaction forces (weak and strong) between particles.

Currently, binary interactions between macro-particles (or super-particles)
are studied under collision schemes that model impact ionization333Process
that generates novel _macro-electrons_ and _macro-ions_ Pukhov (2015). Group-
wise $n$-ary interactions are studied under several _interaction cross-
sections_ Martinez et al. (2019), where each order-$n$ cross section focuses
on a $n$-ary interaction, for a fixed $n$ (e.g., $3$-cross section focuses on
ternary interactions only). Prediction of future interactions in evolving
plasma has garnered a lot of attention in the field of _computational plasma
physics_ Shalaby et al. (2017). The outcome of each near-future prediction of
higher-order $n$-ary interactions can potentially relate to the detection of
novel particles (via novel energy readings on detectors). Additionally, it has
an added advantage of being able to optimally place fewer detectors, which
also results in cost-savings.

###  A.2. Structural biology and biochemistry

Proteins are composed of amino acids linked by covalent peptide bonds.
_Protein structure networks_ (PSNs) Greene (2012) model protein structure and
its mutations as an evolving graph, where an amino acid is considered as a
node and an interaction between a pair of amino acids is a link/edge in the
PSN. It is well-known that the protein structure is directly correlated to the
function of the protein Mannige (2014). Key to the study of protein structure
is gaining an understanding of its structural stability and dynamics. The
higher-order interplay between select _amino acid groups_ within PSNs, modeled
as cliques, are referred to as “protein sectors” Halabi et al. (2009).
Predicting the formation of such sectors in PSNs gives deeper insights into
the robustness of protein structures to mutations. For example, in _protein
therapeutics_ , understanding the protein stability allows them to decide the
shelf-life of a protein drug and their rates of decay (half-lives) in a
patient’s bloodstream.

##  Appendix B Example of Simplex and Related Notions

######  Example 1.

We begin by computing the $k$-balls centered at $1$-simplex $[9,10]$ in
$\mathcal{G}_{t}$ and $\mathcal{G}_{t-1}$, respectively.

The $k$-ball at time $t$ for $k=1$ (i.e., 1-hop vertices only) centered at
$[9,10]$ is:

$B_{t,1}([9,10])=B_{t,1}([9])\cup B_{t,1}([10])$. This is the union of
$k$-balls at underlying vertices $9$ and $10$ according to Definition 3.

$\displaystyle B_{t,1}([9,10])$ $\displaystyle=B_{t,1}([9])\cup B_{t,1}([10])$
$\displaystyle=\\{9,13,8,5,6,10\\}\cup\\{10,14,13,9,6,7,11,15\\}$
$\displaystyle=\\{9,13,8,5,6,10,14,7,11,15\\}$

Similarly, The $k$-ball at previous time step $t-1$ for $k=1$ (i.e., 1-hop
vertices only) centered at $[9,10]$ is:

$\displaystyle B_{t-1,1}([9,10])$ $\displaystyle=B_{t-1,1}([9])\cup
B_{t-1,1}([10])$
$\displaystyle=\\{9,13,8,5,6,10\\}\cup\\{10,14,13,9,6,11,15\\}$
$\displaystyle=\\{9,13,8,5,6,10,14,11,15\\}$

Notice that there is only a difference of vertex $7$ missing from set
$B_{t-1,1}([9,10])$ as compared to set $B_{t,1}([9,10])$ at time $t$.

Now, we calculate the subcomplex spanned by $B_{t,1}([9,10])$ Then444all
simplices are placed on a line each in increasing order of their dimension,

$\displaystyle\mathcal{G}^{\prime}_{t}([9,10])=$
$\displaystyle\\{[9],[13],[8],[5],[6],[10],[14],[11],[15],[7],$ (8)
$\displaystyle[9,10],[9,6],[9,5],[9,8],[9,13],[10,6],[10,7],[10,11],[10,15],[10,14],[10,13],$
$\displaystyle[9,6,10],[10,14,15]\\}$

Now, we compute the compressed $f$-vector notation of
$\mathcal{G}^{\prime}_{t}([9,10])$ to get

$f(\mathcal{G}^{\prime}_{t}([9,10]))=(1,10,11,2)$

Finally, the neighborhood $N_{t}([9,10])=(1,10,11,2)$.

##  Appendix C Further details of the experiment

###  C.1. Datasets

We report results on real-world dynamic graph datasets sourced from Benson et.
al. Benson et al. (2018). Each dataset is a set of _timestamped simplices_
(represented as a set of nodes). In each dataset, let $n$, $m$, and $x$ denote
the number of nodes, edges formed, and timestamped simplices, respectively.
_Enron_ ($n=143$, $m=1.8$K, $x=5$K) and _EU_ ($n=998$, $m=29.3$K, $x=8$K)
model email networks where nodes are email addresses and all recipients of an
email form a simplex in the network. _Contact_ ($n=327$, $m=5.8$K, $x=10$K) is
a proximity graph where nodes represent persons and a simplex is a set of
persons in close proximity to each other. _NDC_ ($n=1.1$K, $m=6.2$K, $x=12$K)
is a drug network from the _National Drug Code_ directory, where nodes are
class labels and a simplex is formed when a set of class labels appear
together on a single drug.

###  C.2. Compared methods

_Adamic-Adar_ Adamic and Adar (2001) and the _Jaccard Coefficient_ Salton and
McGill (1986) measure link probability between two nodes based on the
closeness of their respective feature vectors. _Preferential attachment_
Mitzenmacher (2004) has received considerable attention as a model of growth
of networks as they model future link probability as the product of the
current number of neighbors of the two nodes. Motivated by resource allocation
in transportation networks (much alike the Optimal Transport (OT) problem),
_Resource allocation index_ Zhou et al. (2009b) proposes a node $x$ tries to
transmit a unit resource to node $y$ via common neighbors that play the role
of _transmitters_ and similarity is measured by the amount of the resource $y$
received from $x$. Node2vec Grover and Leskovec (2016b) and SEAL Zhang and
Chen (2018b) are deep-learning based graph embedding methods that are used in
link prediction.

######  Remark 1 (Difference between our setting and simplical closure).

The closest work Benson et al. (2018) proposed predicting a “simplicial
closure” event where at time $t$ there exists a set of nodes which are
pairwise edge connected and the task is to predict whether at time $t+1$ there
will arrive a simplex which covers all these nodes. This phenomenon was termed
as _simplicial closure_. For example, authors $A$, $B$ and $C$ have all co-
authored in pairs (i.e., $\\{A,B\\},\\{A,C\\}$ and $\\{B,C\\}$) and a
simplicial closure event would take place at $t+1$, if a simplex $\\{A,B,C\\}$
arrives, implying that all three authors co-author on a single paper. Our
prediction task significantly diverges and aims to solve a different problem.
Considering our previous example, we are given a single co-authorship
relationship between say $A$ and $B$ at time $t$, we predict whether authors
$A$ and $B$ will co-author with a third author $C$ (ternary co-authorship
relationship) on a single paper at time $t+1$ in the future.

##  Appendix D Proof for Consistency

For preparation, we rewrite the estimator $\widehat{g}_{T}$. Plugging in the
definition of our kernel (Equation 5) into the equation of our estimator
(Equation 6) along with the definitions of $P_{t}(\cdot,\cdot)$ and
$P_{t}^{\tau}(\cdot,\cdot)$ to replace the indicator variables with actual
counts, we obtain the following simplification of Equation 6. Then, it is
reformulated as

$\displaystyle\widetilde{g}_{T}\left(F\right)=\frac{\sum_{t^{\prime}=T-p}^{T}\sum_{\sigma^{(d)}_{j}\in\mathcal{G}_{t^{\prime}}{(d)}}\left(\lvert
P^{\tau}_{t^{\prime}+1}(\sigma^{(d)}_{j},F)\rvert+\beta\sum_{s\in\Gamma(F,\delta)}\lvert
P^{\tau}_{t^{\prime}+1}(\sigma^{(d)}_{j},s)\rvert\right)}{\sum_{t^{\prime}=T-p}^{T}\sum_{\sigma^{(d)}_{j}\in\mathcal{G}_{t^{\prime}}{(d)}}\left(\lvert
P_{t^{\prime}+1}(\sigma^{(d)}_{j},F)\rvert+\beta\sum_{s\in\Gamma(F,\delta)}\lvert
P_{t^{\prime}+1}(\sigma^{(d)}_{j},s)\rvert\right)}.$ (9)

When we set $\beta=0$, we look for other pairs whose feature corresponds to
$F$ and we calculate a fraction of how many such close by pairs actually form
a $(d+1)$-simplex at time $t^{\prime}+1$. This fraction is summed across
various $d$-simplices $\sigma_{j}^{(d)}$ by varying $j$ and also across
various discrete time steps by varying $t^{\prime}$. Setting $\beta>0$ allows
our estimator to _smooth_ over close by features.

It turns out simpler to study a proxy estimator $\widehat{g}_{T}$, which omits
the smoothing of the original one. We will show that it is asymptotically
equivalent to $\widetilde{g}$. Let $\lvert\mathcal{G}_{t}{(d)}\rvert$ denote
the total number of $d$-simplices in GSC $\mathcal{G}_{t}$ at time $t$. Then,
for a feature $F$, we define the proxy estimator as

$\displaystyle\widehat{g}_{T}(F):=\frac{\widehat{h}_{T}(F)}{\widehat{d}_{T}(F)},$
(10)

where the terms are defined as

$\displaystyle\widehat{h}_{T}(F)$
$\displaystyle=\frac{1}{T-p}\sum_{t=p}^{T-1}\sum_{j=1}^{\lvert\mathcal{G}_{t}(d)\rvert}\frac{\lvert
P^{\tau}_{t+1}(\sigma^{(d)}_{j},F)\rvert}{\lvert\mathcal{G}_{t}(d)\rvert},\mbox{~{}and~{}}\widehat{d}_{T}(F)=\frac{1}{T-p}\sum_{t=p}^{T-1}\sum_{j=1}^{\lvert\mathcal{G}_{t}(d)\rvert}\frac{\lvert
P_{t+1}(\sigma^{(d)}_{j},F)\rvert}{\lvert\mathcal{G}_{t}(d)\rvert}.$

Observe that $\sum_{i=1}^{\lvert\mathcal{G}_{t}(d)\rvert}$ ranges over the
total number of $d$-simplices in $\mathcal{G}_{t}$555In practice, the total
number of $d$-simplices is much less than the maximum possible cliques with
$n$ vertices. and note that $\lvert\mathcal{G}_{t}(d)\rvert$ changes with time
step $t$. Recall that the terms $P^{\tau}_{t+1}(\cdot,\cdot)$ and
$P_{t+1}(\cdot,\cdot)$ count all _actual_ and _possible_ formation of
$[\sigma_{i}^{(d)},\widetilde{v}_{m}]$ ($(d+1)$-simplex) that have the _same
feature vector_ $F$. Lemma 1 in the supplementary material proves
$\lvert\widetilde{g}_{T}(F)-\widehat{g}_{T}(F)\rvert\to 0$ as $\beta\to 0$.
First of all, we show the validity of the proxy estimator $\widehat{g}_{T}$.

######  Lemma 1 (Approximation by proxy).

We obtain

$\lvert\widetilde{g}_{T}(F)-\widehat{g}_{T}(F)\rvert=O(\beta),\forall F.$

###### Proof of Lemma 1.

Recall that $\Gamma(F,\delta)$ denotes the set of features at $L_{1}$-distance
at most $\delta$ from $F$. We denote by $\lvert\Gamma(F,\delta)\rvert$ the
cardinality of this set.

$\widetilde{g}_{T}(F)=\frac{\widehat{h}_{T}(F)+C_{1}}{\widehat{d}_{T}(F)+C_{2}}$

where $C_{1}=\beta\sum_{j,t}\sum_{s\in\Gamma(F,\delta)}\lvert
P^{\tau}_{t+1}(\sigma^{(d)}_{j},s)\rvert$ and
$C_{2}=\beta\sum_{j,t}\sum_{s\in\Gamma(F,\delta)}\lvert
P_{t+1}(\sigma^{(d)}_{j},s)\rvert$. Due to the finiteness of features in
$P^{\tau}_{t+1}(\cdot,\cdot)$ and $P_{t+1}(\cdot,\cdot)$, we have that
$C_{1}=C_{2}=O(\beta)$. Both $C_{1}$ and $C_{2}$ are non-negative integers.
So,

$\lvert\widetilde{g}_{T}(F)-\widehat{g}_{T}(F)\rvert=\left\lvert\frac{\widehat{h}_{T}(F)+C_{1}}{\widehat{d}_{T}(F)+C_{2}}-\frac{\widehat{h}_{T}(F)}{\widehat{d}_{T}(F)}\right\rvert=O(\beta)$

In the last step, the second fraction is a positive constant and can thus be
ignored from our asymptotic analysis because both $\widehat{h}_{T}$ and
$\widehat{d}_{T}$ are bounded. ∎

Next, we prove the convergence of the proxy estimator $\widehat{g}_{T}(F)$. As
the first step, we describe the detail of its decomposition in (7). To
simplify the notation, we will drop $F$ from all estimator notations.

######  Proposition 2.

As written in Equation 7, we obtain

$\displaystyle\widehat{g}_{T}-g=\mathcal{V}_{T}+\mathcal{B}_{T}.$

Furthermore, with the event $S_{C}$, there exist a stochastic terms $q_{t}$
for $t$ such as

$\displaystyle\mathcal{V}_{T}=\frac{(T-p)^{-1}\sum_{t=p}^{T-1}q_{t}}{\widehat{d}_{T}}.$

###### Proof of Proposition 2.

With the definition of $B_{T}(F,C)$, we have

$\displaystyle\widehat{g}_{T}-g$
$\displaystyle=\frac{\widehat{h}_{T}}{\widehat{d}_{T}}-g$ (11)
$\displaystyle=\frac{\widehat{h}_{T}-g\widehat{d}_{T}}{\widehat{d}_{T}}$
$\displaystyle=\frac{[\widehat{h}_{T}-g\widehat{d}_{T}]-\mathds{E}[\widehat{h}_{T}-g\widehat{d}_{T}\mid
S_{C}]+\mathds{E}[\widehat{h}_{T}-g\widehat{d}_{T}\mid
S_{C}]}{\widehat{d}_{T}}$
$\displaystyle=\frac{[\widehat{h}_{T}-g\widehat{d}_{T}]-\mathds{E}[\widehat{h}_{T}-g\widehat{d}_{T}\mid
S_{C}]}{\widehat{d}_{T}}+\frac{\mathds{E}[\widehat{h}_{T}\mid
S_{C}]-g\mathds{E}[\widehat{d}_{T}\mid S_{C}]}{\widehat{d}_{T}}$
$\displaystyle={\frac{[\widehat{h}_{T}-g\widehat{d}_{T}]-\mathds{E}[\widehat{h}_{T}-g\widehat{d}_{T}\mid
S_{C}]}{\widehat{d}_{T}}}+\frac{B_{T}(F,C)\mathds{E}[\widehat{d}_{T}\mid
S_{C}]}{\widehat{d}_{T}}$ $\displaystyle=\mathcal{V}_{T}+\mathcal{B}_{T}.$
(12)

We are interested in the asymptotic behavior of the Markov chain at time
$T\rightarrow\infty$. Let $F$ denote
$F_{T}(\sigma_{i}^{(d)},\widetilde{v}_{m})$. Recall that the terms $\lvert
P^{\tau}_{T+1}(\cdot,\cdot)\rvert$ and $\lvert P_{T+1}(\cdot,\cdot)\rvert$
count all _actual_ and _possible_ formation of
$[\sigma_{i}^{(d)},\widetilde{v}_{m}]$ ($(d+1)$-simplex) that result in the
_same feature vector_ $F$ at time $T+1$. Our Markov chain has a _finite state
space_ and hence belongs to a closed communication class with probability
approaching to $1$. We provide a statistical consistency conditional on
$S_{C}$ for any communication class $C$.

For a given time step $t$, we define

$\displaystyle\widehat{h}_{T}(t):=\frac{1}{\lvert\mathcal{G}_{t}(d)\rvert}\sum_{j=1}^{\lvert\mathcal{G}_{t}(d)\rvert}\lvert
P^{\tau}_{t+1}(\sigma_{j}^{(d)},F)\rvert$ (13)
$\displaystyle\widehat{d}_{T}(t):=\frac{1}{\lvert\mathcal{G}_{t}(d)\rvert}\sum_{j=1}^{\lvert\mathcal{G}_{t}(d)\rvert}\lvert
P_{t+1}(\sigma_{j}^{(d)},F)\rvert$

Note that $\widehat{h}_{T}=\frac{1}{T-p}\sum_{t=p}^{T-1}\widehat{h}_{T}(t)$
and $\widehat{d}_{T}=\frac{1}{T-p}\sum_{t=p}^{T-1}\widehat{d}_{T}(t)$.

Let us set

$\displaystyle
q_{t}:=[\widehat{h}_{T}(t)-g\widehat{d}_{T}(t)]-\mathds{E}[\widehat{h}_{T}(t)-g\widehat{d}_{T}(t)\mid
S_{C}]$ (14)

Note that $q_{t}$ is the numerator of the _stochastic term_ in Equation 7 and
a _bounded deterministic function_ of $S_{C}$ at a given time step $t$. ∎

For the stochastic term $\widehat{d}_{T}$ which appears in the denominator, we
show its convergence. The following two lemmas provide the result.

######  Lemma 2.

If the GSC process is $\alpha$-mixing, then, as $T\to\infty$, we obtain

$\displaystyle\operatorname{Var}(\widehat{h}_{T}(F)|S_{C})\to
0,\mbox{~{}and~{}}\operatorname{Var}(\widehat{d}_{T}(F)|S_{C})\to 0,$

for any $F$.

###### Proof of Lemma 2.

We show that variance divided by $T$ converges to a non-negative constant. Let
$U_{T}:=\sum_{t}q_{t}/\sqrt{T}$, where $q_{t}$ (as shown in Equation 14) is a
bounded deterministic function of the state of $X_{t}$ at time $t$. As
demonstrated in Sarkar et. al. Sarkar et al. (2014a), we too break our
weighted sum $U_{T}$ across three time intervals: (i) $[1,T_{C}-1]$, (ii)
$[T_{C},T_{C}+M-1]$, and (iii) $[T_{C}+M,T]$, where $M$ is a constant. Now,
from Sarkar et al. (2014a), we simply apply Lemma $5.7$ to get that
$\mathds{E}[\operatorname{Var}(U_{T}\mid\mathcal{E}(T_{C}),S_{C})\mid
S_{C}]\rightarrow\sigma_{c}$, for some $\sigma_{c}\geq 0$ and from Lemma
$5.8$, we have that
$\operatorname{Var}(\mathds{E}[U_{T}\mid\mathcal{E}(T_{C}),S_{C}]\mid
S_{C})=o(1)$.

Now, since the law of total variance provides

$\operatorname{Var}(U_{T}\mid
S_{C})=\mathds{E}[\operatorname{Var}(U_{T}\mid\mathcal{E}(T_{C}),S_{C})\mid
S_{C}]+\operatorname{Var}(\mathds{E}[U_{T}\mid\mathcal{E}(T_{C}),S_{C}]\mid
S_{C})$

we use the previous results from Lemmas $5.7$ and $5.8$ in Sarkar et al.
(2014a) to get

$\operatorname{Var}(U_{T}\mid S_{C})\rightarrow\sigma_{c}\text{ as
}T\rightarrow 0\text{, for some constant }\sigma_{c}\geq 0$

Plugging in the definition of $q_{t}$ into $U_{T}$ and calculating
$\operatorname{Var}(U_{T}\mid S_{C})$, it follows trivially that
$\operatorname{Var}(\widehat{h}_{T}\mid S_{C})\rightarrow 0$ and
$\operatorname{Var}(\widehat{d}_{T}\mid S_{C})\rightarrow 0$ as
$T\rightarrow\infty$. We refer readers to Remark $5.10$ in Sarkar et al.
(2014a) to see how these results also hold in the case when $C$ is
_aperiodic_. ∎

######  Lemma 3.

If the GSC process is $\alpha$-mixing, then there exist a function $R(C)$ with
a deterministic function of class $C$ denote, such as

$\lim_{T\rightarrow\infty}\mathds{E}[\widehat{d}_{T}(F)\mid\mathcal{E}(T_{C}),S_{C}]=R(C),\text{
and }\lim_{T\rightarrow\infty}\mathds{E}[\widehat{d}_{T}(F)\mid S_{C}]=R(C).$

###### Proof of Lemma 3.

We know by definition that

$\displaystyle\mathds{E}[\widehat{d}_{T}(F)\mid\mathcal{E}(T_{C}),S_{C}]=\frac{1}{T-p}\sum_{t=p}^{T-1}\sum_{j=1}^{|\mathcal{G}_{t}(d)|}\mathds{E}\left[\frac{|P_{t+1}(\sigma^{(d)}_{j},F)|}{|\mathcal{G}_{t}(d)|}\biggm{|}\mathcal{E}(T_{C}),S_{C}\right]$
(15)

This is an average of terms
$\mathds{E}\left[\frac{|P_{t+1}(\sigma^{(d)}_{j},F)|}{|\mathcal{G}_{t}(d)|}\biggm{|}\mathcal{E}(T_{C}),S_{C}\right]$
spanning across $d$-simplices with indices
$j\in\\{1,\cdots,|\mathcal{G}_{t}(d)|\\}$ and discrete time steps
$t\in\\{p,\cdots,T-1\\}$.

For ease of notation, let

$X_{j}:=\frac{|P_{t+1}(\sigma^{(d)}_{j},F)|}{|\mathcal{G}_{t}(d)|}$

$X_{j}$ denotes the total number of possible $(d+1)$-simplices with a $d$-face
as $\sigma^{(d)}_{j}$ divided by the total number of $d$-simplices in
$\mathcal{G}_{t}$.

In the R.H.S. of Equation 15, the term inside the summation is simplified as

$\mathds{E}[X_{j}\mid\mathcal{E}(T_{C}),S_{C}]=\sum_{x}x\mbox{Pr}[X_{j}=x\mid\mathcal{E}(T_{C}),S_{C}]$

We know that both $P_{t+1}(\sigma^{(d)}_{j},F)$ and $\mathcal{G}_{t}(d)$ are
fully determined given the current state $S_{t}$ of the Markov chain. Let
$\mathds{I}_{S}(Y)$ denote an indicator variable of whether “$Y$ is in state
$S$” or not.

We have,

$\mbox{Pr}[X_{j}=x\mid\mathcal{E}(T_{C}),S_{C}]=\sum_{S}\mathds{I}_{S}(X_{j}=x)\mbox{Pr}[S_{t}=S\mid\mathcal{E}(T_{C}),S_{C}]$

As a result, the R.H.S. of Equation 15 becomes

$\displaystyle\frac{1}{T}\sum_{t}\sum_{S}(\sum_{j,x}x\mathds{I}_{S}(X_{j}=x))\mbox{Pr}[S_{t}=S\mid\mathcal{E}(T_{C}),S_{C}]$
(16)

Let $\lambda(S)=\sum_{j,x}x\mathds{I}_{S}(X_{j}=x)$ as this term is fully
determined by state $S$. Then, Equation 16 can be rewritten as

$\displaystyle\sum_{S}\lambda(S)\frac{\sum_{t}\mbox{Pr}[S_{t}=S\mid\mathcal{E}(T_{C}),S_{C}]}{T}$
(17)

Due to stationarity, the average
$\sum_{t}\mbox{Pr}[S_{t}=S\mid\mathcal{E}(T_{C}),S_{C}]/T$ will converge to a
constant function of state $S$, denoted by $R(S)$. Given that $\lambda(S)$ is
bounded and the average term converges to a constant $R(S)$, we say that
Equation 17 converges to some constant $R(C)>0$, where $R(C)$ is a
deterministic function of communication class $C$. This proves the first part.

A simple application of the _tower property of expectation_ followed by the
_dominated convergence theorem_ shows that
$\lim_{T\rightarrow\infty}\mathds{E}[\widehat{d}(F)\mid S_{C}]=R_{C}$. This
completes the proof for the result. ∎

Then, we are ready to prove the convergence of the variance term.

######  Proposition 3 (Variance).

If the GSC filtration process is $\alpha$-mixing, then, conditional on
$S_{C}$, we obtain $\mathcal{V}_{T}\overset{p}{\to}0$ as $T\to\infty$.

###### Proof of Proposition 3.

By Proposition 2, the term $\mathcal{V}_{T}$ is written as
$\frac{(T-p)^{-1}\sum_{t}q_{t}}{\widehat{d}_{T}}$.

For the denominator, Lemma 3 shows that $\mathds{E}[\widehat{d}_{T}\mid
S_{C}]\rightarrow R(C)$, where $R(C)$ is a positive deterministic function of
class $C$. Also, Lemma 2 states that $Var(\widehat{d}_{T}\mid
S_{C})\rightarrow 0$ holds as $T\rightarrow\infty$. Thus,
$\widehat{d}_{T}\overset{p}{\to}R(C)>0$ holds. Also, $\mathcal{V}_{T}$ is
asymptotically well defined for class $C$.

For the nominator, Lemma 2 also shows that
$\lim_{T\rightarrow\infty}Var(\sum_{t}q_{t}/T\mid S_{C})=0$ as
$T\rightarrow\infty$ and $\mathds{E}[q_{t}\mid S_{C}]=0$, therefore we have
$\lim_{T\rightarrow\infty}1/T\sum_{t}q_{t}\overset{qm}{\to}0$ conditioned on
$S_{C}$.

By the continuous mapping theorem, we obtain the statement. ∎

Next, we discuss the bias term $\mathcal{B}_{T}$. To this aim, we rewrite the
term $B_{T}(F,C)$ as follows:

$B_{T}(F,C)=\frac{(\mathds{E}[\widehat{h}_{T}(F)\mid
S_{C}]-g(F)\mathds{E}[\widehat{d}_{T}(F)\mid
S_{C}])}{\mathds{E}[\widehat{d}_{T}(F)\mid S_{C}]}.$

######  Lemma 4.

Given that Assumption 1 holds and that when $T\rightarrow\infty$, the
bandwidth parameter $b\rightarrow 0$. Then, we have that
$B_{T}(F,C)=O(b)=o(1)$ as $T\to\infty$.

###### Proof of Lemma 4.

For $t\in[p,T-1]$, $j\in[1,|\mathcal{G}_{t}(d)|]$ and a fixed feature vector
$F$, the numerator of $B_{T}(F,C)$ can be expressed as an average of terms of
the form

$\displaystyle
A_{t}:=\mathds{E}\left[\frac{|P^{\tau}_{t+1}(\sigma^{(d)}_{j},F)|}{|\mathcal{G}_{t}(d)|}\biggm{|}S_{C}\right]-\mathds{E}\left[\frac{|P_{t+1}(\sigma^{(d)}_{j},F)|}{|\mathcal{G}_{t}(d)|}\biggm{|}S_{C}\right]g(F)$
(18)

The first term in Equation 18 can be rewritten using the _tower property_ as

$\mathds{E}\left[\mathds{E}\left[\frac{|P^{\tau}_{t+1}(\sigma^{(d)}_{j},F)|}{|\mathcal{G}_{t}(d)|}\biggm{|}\mathcal{E}(T_{C}),S_{C}\right]\biggm{|}S_{C}\right]$

When we condition on $\mathcal{E}(T_{C})$, it makes
$\frac{|P^{\tau}_{t+1}(\sigma^{(d)}_{j},F)|}{|\mathcal{G}_{t}(d)|}$
conditionally independent of $S_{C}$, if $t>T_{C}$. Also, for $t\geq T_{C}$,
we have

$\displaystyle\mathds{E}\left[\frac{|P^{\tau}_{t+1}(\sigma^{(d)}_{j},F)|}{|\mathcal{G}_{t}(d)|}\biggm{|}\mathcal{E}(T_{C}),S_{C}\right]=\frac{|P_{t+1}(\sigma^{(d)}_{j},F)|}{|\mathcal{G}_{t}(d)|}g(F_{t}(\sigma^{(d)}_{j},\widetilde{v}_{n}))$
(19)

where $\widetilde{v}_{n}\in B_{t-1,k}(\sigma^{(d)}_{j})$.

Given the result in Equation 19 and the fact that the term
$\frac{|P^{\tau}_{t+1}(\sigma^{(d)}_{j},F)|}{|\mathcal{G}_{t}(d)|}$ is bounded
results in

$\displaystyle\mathds{E}\left[\frac{|P^{\tau}_{t+1}(\sigma^{(d)}_{j},F)|}{|\mathcal{G}_{t}(d)|}\biggm{|}\mathcal{E}(T_{C}),S_{C}\right]$
(20)
$\displaystyle\leq\frac{|P_{t+1}(\sigma^{(d)}_{j},F)|}{|\mathcal{G}_{t}(d)|}g(F_{t}(\sigma^{(d)}_{j},\widetilde{v}_{n}))\mathds{I}[T_{C}\leq
t]+c\mathds{I}[T_{C}>t]$
$\displaystyle\leq\frac{|P_{t+1}(\sigma^{(d)}_{j},F)|}{|\mathcal{G}_{t}(d)|}g(F_{t}(\sigma^{(d)}_{j},\widetilde{v}_{n}))+c\mathds{I}[T_{C}>t],$

where $c>0$ an existing constant.

Now, the numerator of $B_{T}(F,C)$ can be upper bounded as

$\displaystyle\left\lvert\sum_{t}A_{t}/T\right\rvert\leq$
$\displaystyle\sum_{t}\frac{1}{T}\left\lvert\mathds{E}\left[\frac{|P_{t+1}(\sigma^{(d)}_{j},F)|}{|\mathcal{G}_{t}(d)|}(g(F_{t}(\sigma^{(d)}_{j},\widetilde{v}_{n}))-g(F))\biggm{|}S_{C}\right]\right\rvert$
(21) $\displaystyle+c^{\prime}\sum_{t}\mbox{Pr}[T_{C}>t]/T.$

The second term in Equation 21 vanishes as $T\rightarrow\infty$ because it is
of order $O(\mathds{E}[T_{C}]/T)$. Thus, the numerator of $B_{T}(F,C)$ is an
average of terms of the form

$\displaystyle\mathds{E}\left[\frac{|P_{t+1}(\sigma^{(d)}_{j},F)|}{|\mathcal{G}_{t}(d)|}(g(F_{t}(\sigma^{(d)}_{j},\widetilde{v}_{n}))-g(F))\biggm{|}S_{C}\right].$
(22)

Our feature vector counts and simplex neighborhoods are finite because
$|\mathcal{G}_{t}(d)|$ is bounded. The expectation in Equation 22 is just a
summation of finite terms. We set
$F^{\prime}=F_{t}(\sigma^{(d)}_{j},\widetilde{v}_{n})$, and make use of our
smoothness assumption 1, so that

$\lvert g(F^{\prime})-g(F)\rvert=O(\kappa(-\|F-F^{\prime}\|_{1}/b)).$

We also use Lemma 3 to say that the denominator of our bias term converges to
a constant $R(C)$. So,

$\displaystyle B_{T}(F,C)$
$\displaystyle=O(\kappa(-\|F-F^{\prime}\|_{1}/b))=O(b).$

The last equality follows the property of $\kappa$ in the Schwartz space. Now,
since $B_{T}(F,C)=O(b)$ and $b\rightarrow 0$ as $T\rightarrow\infty$, then we
have that $B_{T}(F,C)=o(1)$. This completes the proof. ∎

Then, we prove the convergence of the bias. Then, we have the result:

######  Proposition 4 (Bias).

If Assumption 1 holds, then, conditional on $S_{C}$,
$\mathcal{B}_{T}\overset{p}{\to}0$ holds as $T\to\infty$.

###### Proof of Proposition 4.

Proposition 2 shows that $\mathcal{B}_{T}(F)=B_{T}(F,C)/\widehat{d}_{T}(F)$.

By Lemma 2 and 3, we obtain $\widehat{d}_{T}\to R(C)>0$ as $T\to\infty$, as
similarly shown in the proof of Proposition 3. Further, Lemma 4 states that
$B_{T}(F,C)=o(b)$. Combining the results, we obtain the statement. ∎

Now, we can prove the consistency (Theorem 1).

###### Proof of Theorem 1.

For the result of $\widehat{g}_{T}$, we apply the results of Proposition 3 and
4 to the decomposition in the equation 7, then obtain the statement.

For the result of $\widetilde{g}_{T}$, we additionally combine the result of
Lemma 1, then obtain the statement. ∎

##  Appendix E Proof for Asymptotic Normality

###  E.1. Introduction to Wasserstein distance and approximation technique

We denote by $BL(\mathds{R})$ the space of such bounded functions $h$ that are
$1$-Lipschitz. More formally,

$\lVert h\rVert_{\infty}=\sup_{x\in\mathds{R}}\lvert h(x)\rvert<\infty\text{
and }Lip(h)=1\text{, where }Lip(h)=\sup_{x\neq y}\frac{|h(x)-h(y)|}{|x-y|}$

So, $h\in BL(\mathds{R})$.

We now make use of the Wasserstein metric to measure the distance between
distributions. Therefore, our estimator represented as $W$ can be shown to
converge to $Z$, when the Wasserstein distance between $W$ and $Z$’s
underlying distributions converges to zero. We have that

$\displaystyle d_{w}(W,Z)=\sup_{h\in
BL(\mathds{R})}\lvert\mathds{E}h(W)-\mathds{E}h(Z)\rvert$ (23)

####  E.1.1. Introduction to Stein’s Method for normal approximation

Stein Stein et al. (1986) introduced a powerful technique to estimate the rate
of convergence of sums of weakly dependent r.v.s to the standard normal
distribution. A remarkable feature of Stein’s method is that it can be applied
in many circumstances where dependence plays a role, therefore we propose an
adaptation of Stein’s method to our setting of dynamic GSCs.

Given a standard normal r.v. $Z$, Stein’s lemma (stated below) provides a
characterization of $Z$’s distribution.

######  Lemma 5 (Stein’s Lemma Chen et al. (2010)).

If $W$ has a standard normal distribution, then

$\displaystyle\mathds{E}f^{\prime}(W)=\mathds{E}[Wf(W)],$ (24)

for all absolutely continuous functions $f:\mathds{R}\rightarrow\mathds{R}$
with $\mathds{E}\lvert f^{\prime}(Z)\rvert<\infty$. Conversely, if Equation 24
holds for all bounded, continuous, and piecewise continuously differentiable
functions $f$ with $\mathds{E}\lvert f^{\prime}(Z)\rvert<\infty$, then $W$ has
a standard normal distribution.

In order to show that a r.v. $W$ has a distribution _close to_ that of a
target distribution of $Z$, one must compare the values of expectations of the
two distributions on some collection of bounded functions
$h:\mathds{R}\rightarrow\mathds{R}$. Here, Stein’s lemma (Lemma 5) shows that
$W\stackrel{{\scriptstyle d}}{{=}}Z$, if

$\displaystyle\mathds{E}f^{\prime}(W)-\mathds{E}[Wf(W)]=0$ (25)

holds. Observe that if the distribution of $W$ is close to that of $Z$’s
distribution, then evaluating the L.H.S of Equation 25 when $W$ is replaced by
$Z$ would result in a small value. Putting these difference equations
together, the following linear differential equation known as _Stein’s
equation_ is arrived at

$\displaystyle f^{\prime}(W)-Wf(W)=h(W)-\mathds{E}h(Z)$ (26)

The $f$ that satisfy Equation 26 with $h\in BL(\mathds{R})$ must satisfy the
following conditions for all $y,z\in\mathds{R}$

$\displaystyle\lVert f\rVert\leq 2\text{ , }\lVert f^{\prime}\rVert\leq
2\text{ , }\lVert f^{\prime\prime}\rVert\leq\sqrt{2/\pi}$ (27)
$\displaystyle|f^{\prime}(y+z)-f^{\prime}(y)|\leq D|z|$

where $h_{0}(y)=h(y)-\mathds{E}h(Z)$, $c_{1}=\sup_{x\geq 0}\xi(x)$,
$c_{2}=\sup_{x\geq 0}x(1-x\xi(x))$, and $\xi(x)=(1-\Phi)/\phi$ (where
$\Phi(x)$ is the _distribution function_ and $\phi(x)=\Phi^{\prime}(x)$).
Then, $D=(c_{1}+c_{2})\lVert h_{0}\rVert_{\infty}+2$ is a constant.
Additionally, we have a bound on the covariance, given the dependent r.v.’s
are also bounded.

If $\displaystyle\mbox{Pr}\\{|X|\leq C_{1}\\}=\mbox{Pr}\\{|Y|\leq
C_{2}\\}=1\text{, then}$ (28) $\displaystyle|\operatorname{Cov}(X,Y)|\leq
4C_{1}C_{2}\alpha(r)$

We take a similar approach to Sarkar et. al. Sarkar et al. (2014a) in terms of
using the Wasserstein distance to bound the normal approximation. We first
define the dependency in our GSCs and then propose a notion of _$\alpha$
mixing_ in our context of GSCs. We obtain a tighter bound than the bound
proposed in Sarkar et al. (2014a), by instead following an approach proposed
by Sunklodas Sunklodas (2007).

###  E.2. Gaussian approximation for dependent variables with GSC

In our model, we assume the r.v. $A_{i}$ to represent a $d$-simplex
$\sigma^{(d)}_{i}$ in a GSC $\mathcal{G}$. In order to have a notion of
$\alpha$-mixing in our setting, we must first define a distance between two
$d$-simplices $\sigma_{i}$ and $\sigma_{j}$. We drop the $(d)$ superscript for
brevity and ease of notation. We define this distance as the _Hausdorff_
distance between simplices as

$\displaystyle
d_{H}(\sigma_{i},\sigma_{j})=\max\left\\{\sup_{v\in\sigma_{i}}\inf_{v^{\prime}\in\sigma_{j}}d_{g}(v,v^{\prime}),\sup_{v^{\prime}\in\sigma_{j}}\inf_{v\in\sigma_{i}}d_{g}(v,v^{\prime})\right\\}$
(29)

where $d_{g}(v,v^{\prime})$ counts the number of edges in the geodesic
connecting vertices $v$ and $v^{\prime}$ in $\mathcal{G}^{(1)}_{-}$.

In Stein’s method, the sum of dependent r.v.’s is studied by breaking the sum
$Z_{n}$ into two sets based on the $r$ in mixing coefficient $\alpha(r)$. In
our setting, given a fixed $d$-simplex $\sigma_{i}$, we study two partial sums
pertaining to: 1) all $d$-simplices that are at most $r$-apart from $A_{i}$
and 2) the remaining partial sum after removing the variables pertaining to 1)
from $Z_{n}$.

With this notion of distance between sets of r.v.’s, we modify with slight
deviations from _proposition $4$_ in Sunklodas et. al. Sunklodas (2007) to
accommodate our $\alpha$-mixing in Markov chains based on GSCs. For a sequence
of r.v.’s $X_{1},X_{2},\cdots$ satisfying the $\alpha$-mixing condition, we
write

$Z_{n}=\sum_{i=1}^{n}A_{i}\text{,\hskip 10.00002pt
}A_{i}=\frac{X_{i}}{B_{n}}\text{,\hskip 10.00002pt
}B_{n}^{2}=\mathds{E}(\sum_{i=0}^{n}X_{i})^{2}$

We assume $B_{n}>0$.

$T_{i}^{(m)}$ denotes the contribution of $d$-simplices that are further than
$m$ away from $\sigma_{i}$ and $x(\sigma_{i},r)$ denotes the partial sum of
r.v.’s representing simplices that are exactly $r$ away from $A_{i}$.
Therefore, $\sum_{r=0}^{m}x(\sigma_{i},r)$ gets us all those $d$-simplices
that are greater than or equal to $m$ away from $\sigma_{i}$. We are
interested in the contribution of simplices $r$ away from $\sigma_{i}$ as we
vary $r$ from $0$ to $m$. With Proposition 5, we proceed to derive an upper
bound on $d_{w}(Z_{n},N)$ (i.e., the Wasserstein distance between $Z_{n}$ and
$N$).

######  Proposition 5.

Let $S(\sigma_{i},r)$ denote the set of $d$-simplices whose Hausdorff distance
equals $r$. More formally,

$S(\sigma_{i},r)=\\{\sigma_{j}:{d_{H}(\sigma_{i},\sigma_{j})}=r\\}$

Additionally, let $\widehat{X}$ denote a _mean-centered_ version of r.v. $X$.
Then,

$\displaystyle x(\sigma_{i},r)=\sum_{p\in S(\sigma_{i},r)}A_{p}\text{\hskip
30.00005pt}(x(\sigma_{i},0)=A_{i})$ $\displaystyle
T_{i}^{(m)}=Z_{n}-\sum_{r=0}^{m}x(\sigma_{i},r)\text{,\hskip
10.00002pt}m=0,1,\cdots\text{,\hskip 20.00003pt}(T_{i}^{(-1)}=Z_{n})$

Suppose that $\mathds{E}Z_{n}=0$, $\mathds{E}Z_{n}^{2}=1$, and
$\mathds{E}A_{i}^{2}<\infty$ for all $i=1,\cdots,n$. Let $\epsilon$ be a r.v.
uniformly distributed in $[0,1]$ and independent of other r.v.’s. Let
$f:\mathds{R}\rightarrow\mathds{R}$ be a differentiable function such that
$\sup_{x\in\mathds{R}}|f^{\prime}(x)|<\infty$. Then we have

$\mathds{E}f^{\prime}(Z_{n})-\mathds{E}Z_{n}f(Z_{n})=E_{1}+\cdots+E_{7}$

where

$\displaystyle E_{1}$ $\displaystyle=-\sum_{i=1}^{n}\sum_{r\geq
1}\mathds{E}A_{i}x(\sigma_{i},r)\left[f^{\prime}(T_{i}^{(r)}+\epsilon
x(\sigma_{i},r))-f^{\prime}(T_{i}^{(r)})\right]$ $\displaystyle E_{2}$
$\displaystyle=-\sum_{i=1}^{n}\mathds{E}A_{i}^{2}\left[f^{\prime}(T_{i}^{(0)}+\epsilon
A_{i})-f^{\prime}(T_{i}^{(0)})\right]$ $\displaystyle E_{3}$
$\displaystyle=-\sum_{i=1}^{n}\sum_{r\geq
1}\sum_{q=r+1}^{2r}\mathds{E}\widehat{A_{i}x(\sigma_{i},r)}\delta_{i}^{(q)}\text{
, }E_{4}=-\sum_{i=1}^{n}\sum_{r\geq 1}\sum_{q\geq
2r+1}\mathds{E}\widehat{A_{i}x(\sigma_{i},r)}\delta_{i}^{(q)}$ $\displaystyle
E_{5}$ $\displaystyle=\sum_{i=1}^{n}\sum_{r\geq
1}\mathds{E}A_{i}x(\sigma_{i},r)\sum_{q=0}^{r}\mathds{E}\delta_{i}^{(q)}\text{,
}E_{6}=-\sum_{i=1}^{n}\sum_{q\geq
1}\mathds{E}\widehat{A_{i}^{2}}\delta_{i}^{(q)}\text{,
}E_{7}=\sum_{i=0}^{n}\mathds{E}A_{i}^{2}\mathds{E}\delta_{i}^{(q)},$

and

$\displaystyle\delta_{i}^{(q)}=f^{\prime}(T_{i}^{(q-1)})-f^{\prime}(T_{i}^{(q)}).$

######  Theorem 3.

Consider a sequence of r.v.’s $X_{1},X_{2},\cdots$ that satisfy
$\alpha$-mixing condition (Definition 5). Let $\mathds{E}X_{i}=0$,
$\mbox{Pr}\\{|X_{i}|\leq L\\}=1$, for $i=1,\cdots,n$, for some constant $L>0$.
Then, for every $h\in BL(\mathds{R})$,

$d_{w}(Z_{n},N)\leq
C(S,D)\left(\sum_{i=0}^{n}\mathds{E}|A_{i}|^{3}+\frac{nL^{3}}{B_{n}^{3}}\sum_{r=1}^{n-1}r\alpha(r)\right)$

where $C(S,D)$ is a finite constant which depends on $D$ and
$\max_{i}|S(\sigma_{i},r)|$.

###### Proof of Theorem 3.

Given a $d$-simplex $\sigma_{i}$ and its corresponding r.v. $A_{i}$, recall
that $S(\sigma_{i},r)$ denotes the set of $d$-simplices that are at Hausdorff
distance $r$ away from $\sigma_{i}$. We additionally define $S_{m}$ to denote
the maximum cardinality of $S(\sigma_{i},r)$ for all $i$ and a fixed $r$,
i.e.,

$S_{m}=\max_{i}|S(\sigma_{i},r)|$

In order to upper bound the Wasserstein distance between $Z_{n}$ and $N$, we
estimate the difference $\mathds{E}h(Z_{n})-\mathds{E}h(N)$ using Proposition
5. It was shown in Proposition 5 that this difference is a sum of terms
$E_{1},\cdots,E_{7}$. We will proceed by individually bounding each term.

Bounding $E_{1}$:

$\displaystyle|E_{1}|=\left|\sum_{i=1}^{n}\sum_{r\geq
1}\mathds{E}A_{i}x(\sigma_{i},r)\underbrace{\left[f^{\prime}(T_{i}^{(r)}+\epsilon
x(\sigma_{i},r))-f^{\prime}(T_{i}^{(r)})\right]}_{(i)}\right|$ (30)

We can upper bound term $(i)$ in Equation 30 using Equation 27 by

$D|\epsilon||x(\sigma_{i},r)|\leq D|x(\sigma_{i},r)|\text{ \hskip 20.00003pt
(since $|\epsilon|\leq 1$)}$

Then,

$\displaystyle\left|x(\sigma_{i},r)\left[f^{\prime}(T_{i}^{(r)}+\epsilon
x(\sigma_{i},r))-f^{\prime}(T_{i}^{(r)})\right]\right|\leq
D(x(\sigma_{i},r))^{2}$ (31)

Now, we upper bound $x(\sigma_{i},r)$ as

$\displaystyle x(\sigma_{i},r)\leq S_{m}\left(\frac{L}{B_{n}}\right)$ (32)

since each normalized r.v. in $Z_{n}$ is upper bounded by $L/B_{n}$. The
L.H.S. of Equation 31 is upper bounded by
$DS_{m}^{2}\left(\frac{L^{2}}{B_{n}^{2}}\right)$.

We know that

$\displaystyle\left|\operatorname{Cov}\left(\underbrace{A_{i}},\underbrace{x(\sigma_{i},r)\left[f^{\prime}(T_{i}^{(r)}+\epsilon
x(\sigma_{i},r))-f^{\prime}(T_{i}^{(r)})\right]}\right)\right|$ (33)
$\displaystyle=\underbrace{\mathds{E}\left[A_{i}x(\sigma_{i},r)\left[f^{\prime}(T_{i}^{(r)}+\epsilon
x(\sigma_{i},r))-f^{\prime}(T_{i}^{(r)})\right]\right]}_{(ii)}$
$\displaystyle+\underbrace{\mathds{E}A_{i}}_{=0}\mathds{E}x(\sigma_{i},r)\left[f^{\prime}(T_{i}^{(r)}+\epsilon
x(\sigma_{i},r))-f^{\prime}(T_{i}^{(r)})\right]$

Notice that term $(ii)$ is nothing but the summand in Equation 30. We apply
the covariance bounds (Equation 28), to obtain

$\displaystyle(ii)$ $\displaystyle\leq
4\left(\frac{L}{B_{n}}\right)\left(DS_{m}^{2}\frac{L^{2}}{B_{n}^{2}}\right)\alpha(r)\leq
4DS_{m}^{2}\frac{L^{3}}{B_{n}^{3}}\alpha(r)$ (34)

Then, we have that

$\displaystyle|E_{1}|$ $\displaystyle\leq\sum_{i=1}^{n}\sum_{r\geq
1}\left(4DS_{m}^{2}\frac{L^{3}}{B_{n}^{3}}\alpha(r)\right)\leq
4DS_{m}^{2}\frac{nL^{3}}{B_{n}^{3}}\sum_{r=1}^{n-1}\alpha(r).$ (35)

Here, the summation $\sum_{r=1}^{n-1}$ appears, because there are only $n$
variables that should measure the dependence between each other.

Bounding $E_{2}$:

$\displaystyle|E_{2}|$
$\displaystyle=\left|\sum_{i=1}^{n}\mathds{E}A_{i}^{2}\underbrace{\left[f^{\prime}(T_{i}^{(0)}+\epsilon
A_{i})-f^{\prime}(T_{i}^{(0)})\right]}_{(iii)}\right|$ (36)

Using Equation 27, we have that term $(iii)\leq D|A_{i}|$. Then, $|E_{2}|$ can
simply be bounded as

$\displaystyle|E_{2}|\leq D\sum_{i=1}^{n}\mathds{E}|A_{i}|^{3}$ (37)

Bounding $E_{3}$:

$\displaystyle|E_{3}|$ $\displaystyle=\left|\sum_{i=1}^{n}\sum_{r\geq
1}\sum_{q=r+1}^{2r}\mathds{E}\widehat{A_{i}x(\sigma_{i},r)}\delta_{i}^{(q)}\right|$
(38) $\displaystyle\leq\sum_{i=1}^{n}\sum_{r\geq
1}\sum_{q=r+1}^{2r}\underbrace{\left|\mathds{E}A_{i}x(\sigma_{i},r)\delta_{i}^{(q)}\right|}_{(iv)}+\sum_{i=1}^{n}\sum_{r\geq
1}\sum_{q=r+1}^{2r}\underbrace{\left|\mathds{E}A_{i}x(\sigma_{i},r)\right|\mathds{E}\left|\delta_{i}^{(q)}\right|}_{(v)}$

Let us focus on bounding terms $(iv)$ and $(v)$, separately.

_For term $(iv)$_: We can further split it as

$\displaystyle\left|\mathds{E}\underbrace{A_{i}}_{a}\underbrace{x(\sigma_{i},r)\delta_{i}^{(q)}}_{b}\right|$
(39)

We have previously worked out the bounds for terms $A_{i}$ and
$x(\sigma_{i},r)$. Therefore, we now focus our attention on bounding
$\delta_{i}^{(q)}$.

$\displaystyle\delta_{i}^{(q)}$
$\displaystyle=f^{\prime}(T_{i}^{(q-1)})-f^{\prime}(T_{i}^{(q)})$ (40)
$\displaystyle=f^{\prime}(T_{i}^{(q)}+x(\sigma_{i},r))-f^{\prime}(T_{i}^{(q)})$
$\displaystyle\leq D|x(\sigma_{i},r)|\text{ \hskip 80.00012pt (using
Equation~{}\ref{eq:f_ineq}) }$ $\displaystyle\leq
D\left(S_{m}\frac{L}{B_{n}}\right)$

Therefore, term $b$ (in term $(iv)$) is upper bounded by
$DS_{m}^{2}\frac{L^{2}}{B_{n}^{2}}$.

Applying the covariance bound (Equation 28), we have that

$\displaystyle(iv)$ $\displaystyle\leq
4\left(\frac{L}{B_{n}}\right)\left(DS_{m}^{2}\frac{L^{2}}{B_{n}^{2}}\right)\alpha(r)\leq
4DS_{m}^{2}\frac{L^{3}}{B_{n}^{3}}\alpha(r)$ (41)

_For term $(v)$_: We have calculated some bounds previously, so we can again
split $(v)$ as

$\displaystyle\underbrace{\left|\mathds{E}\underbrace{A_{i}}\underbrace{x(\sigma_{i},r)}\right|}_{c}\mathds{E}\underbrace{\left|\delta_{i}^{(q)}\right|}_{d}$
$\displaystyle\leq\underbrace{\left[4\left(\frac{L}{B_{n}}\right)\left(S_{m}\frac{L}{B_{n}}\right)\alpha(r)\right]}_{\text{using
Eqn~{}\ref{eq:cov} on
$c$}}\underbrace{\left[D\left(S_{m}\frac{L}{B_{n}}\right)\right]}_{\text{using
Eqn~{}\ref{eq:delta} on $d$}}$ (42) $\displaystyle\leq
4DS_{m}^{2}\frac{L^{3}}{B_{n}^{3}}\alpha(r)$

Now, combining the inequalities for terms $(iv)$ and $(v)$, we have that

$\displaystyle|E_{3}|$ $\displaystyle\leq\sum_{i=1}^{n}\sum_{r\geq
1}\sum_{q=r+1}^{2r}8DS_{m}^{2}\frac{L^{3}}{B_{n}^{3}}\alpha(r)\leq
8DS_{m}^{2}\frac{nL^{3}}{B_{n}^{3}}\sum_{r=1}^{n-1}r\alpha(r)$ (43)

Bounding $E_{4}$: We have that

$\displaystyle|E_{4}|$ $\displaystyle=\left|\sum_{i=1}^{n}\sum_{r\geq
1}\sum_{q\geq
2r+1}\mathds{E}\widehat{A_{i}x(\sigma_{i},r)}\delta_{i}^{(q)}\right|$ (44)
$\displaystyle\leq 4DS_{m}^{2}\frac{nL^{3}}{B_{n}^{3}}\sum_{r\geq
1}\sum_{q\geq 2r+1}\alpha(q-r)$ $\displaystyle\leq
4DS_{m}^{2}\frac{nL^{3}}{B_{n}^{3}}\sum_{r=1}^{n-1}r\alpha(r)$

because $\sum_{r\geq 1}\sum_{q\geq
2r+1}\alpha(q-r)<\sum_{r=1}^{n-1}r\alpha(r)$.

Now, using our previous bounds, the remaining terms, i.e., $|E_{5}|$,
$|E_{6}|$, and $|E_{7}|$ are bounded as follows:

$\displaystyle|E_{5}|$ $\displaystyle\leq
4DS_{m}^{2}\frac{nL^{3}}{B_{n}^{3}}\sum_{r=1}^{n-1}r\alpha(r)$ (45)
$\displaystyle|E_{6}|$ $\displaystyle\leq
4DS_{m}^{2}\frac{nL^{3}}{B_{n}^{3}}\sum_{r=1}^{n-1}\alpha(r)$
$\displaystyle|E_{7}|$ $\displaystyle\leq
D\sum_{i=0}^{n}\mathds{E}|A_{i}|^{3}$

Finally, we combine these bounds to arrive at our final upper bound as

$\displaystyle|\mathds{E}f^{\prime}(Z_{n})-\mathds{E}Z_{n}f(Z_{n})|\leq
C(S_{m},D)\left(\sum_{i=0}^{n}\mathds{E}|A_{i}|^{3}+\frac{nL^{3}}{B_{n}^{3}}\sum_{r=1}^{n-1}r\alpha(r)\right)$
(46)

where $C(S_{m},D)$ is a constant term depending on constants $K_{m}$ and $D$.
This completes our proof. ∎

###  E.3. Deferred Proof

In this section, we establish our estimator’s result on the asymptotic
normality.

Let $n=|\mathcal{G}_{t}(d)|$ denote the total number of $d$-simplices in GSC
$\mathcal{G}_{t}$ at time $t$. Recall from Equation 13 that

$\displaystyle\widehat{h}_{T}(t):=\frac{1}{n}\sum_{j=1}^{n}\lvert
P^{\tau}_{t+1}(\sigma_{j}^{(d)},F)\rvert$
$\displaystyle\widehat{d}_{T}(t):=\frac{1}{n}\sum_{j=1}^{n}\lvert
P_{t+1}(\sigma_{j}^{(d)},F)\rvert$

Also, from Equation 14, we have

$q_{t}:=[\widehat{h}_{T}(t)-g\widehat{d}_{T}(t)]-\mathds{E}[\widehat{h}_{T}(t)-g\widehat{d}_{T}(t)\mid
S_{C}]$

Additionally, let us define a variable $p_{t}$ for convenience as follows

$p_{t}:=[\widehat{h}_{T}(t)-g\widehat{d}_{T}(t)]-\mathds{E}[\widehat{h}_{T}(t)-g\widehat{d}_{T}(t)\mid\mathcal{E}(T_{C}),S_{C}]$

Note that $p_{t}$ is conditioned on both $\mathcal{E}(T_{C})$ and $S_{C}$,
while $q_{t}$ is just conditioned on $S_{C}$. Keeping these expressions in
mind, we begin by showing the weak convergence with $p_{t}$.

######  Lemma 6.

Under Assumption 1, given $p_{t}$ for any finite $T_{C}$ and $\sigma_{c}>0$

$\sum_{t\geq
T_{C}+M}p_{t}/\sqrt{T}\overset{d}{\to}\mathcal{N}(0,\sigma_{c}^{2})\text{\hskip
10.00002pt conditioned on }\mathcal{E}(T_{C})\cap S_{C}\text{ as
}T\rightarrow\infty$

###### Proof of Lemma 6.

Given a normalized r.v. $W_{T}$ which is a sum of weakly dependent r.v.’s as

$W_{T}:=\frac{\sum_{t\geq
T_{C}+M}p_{t}}{\sqrt{\operatorname{Var}\left(\sum_{t\geq
T_{C}+M}p_{t}\mid\mathcal{E}(T_{C}),S_{C}\right)}}$

where $p_{t}$ is already _bounded_ and _mean-centered_. We have to show that
our upper bound in Theorem 3 converges to zero, so that according to Stein’s
lemma 5, we have

$W_{T}\overset{d}{\to}\mathcal{N}(0,1)\text{\hskip 10.00002pt conditioned on
}\mathcal{E}(T_{C})\cap S_{C}$

Recall that by Lemma $5.7$ in Sarkar et al. (2014a), we have

$\operatorname{Var}\left(\sum_{t\geq
T_{C}+M}p_{t}\mid\mathcal{E}(T_{C}),S_{C}\right)/T\rightarrow\sigma_{c}^{2}$

Now, we show that conditioned on event $\mathcal{E}(T_{C})\cap S_{C}$, our
bound in Theorem 3 has a convergence rate of $O(T^{-1/2})$.

Note that our $p_{t}$ corresponds to $A_{i}$ and $T$ to $n$ in Theorem 3. As
$p_{t}$ is a function of $S_{t}$, it involves $p+1$ GSCs
$(\mathcal{G}_{t-p+1},\cdots,\mathcal{G}_{t+1})$ and the distance between the
$i$ and $j$-th GSC is defined as $dist(i,j)=\max(|i-j|-(p+1),0)$. Therefore,
you will observe that we have for $x(\sigma_{i},r)$ only $2$ states that are
at distance $r$ apart from $\sigma_{i}$, i.e., $A_{i}-r$ and $A_{i}+r$.
Therefore, $S_{m}=O(1)$.

Given that $\mbox{Pr}\\{|p_{t}|\leq L\\}=1$, for $t=1,\cdots,T$ and for some
constant bound $L>0$, we have

$\displaystyle\sum_{t=1}^{T}\mathds{E}|p_{t}|^{3}+T\frac{L^{3}}{B_{T}^{3}}\sum_{r=1}^{T-1}r\alpha(r)\text{\hskip
20.00003pt (ignoring some constant terms) }$ (47) $\displaystyle\leq
T\frac{L^{3}}{B_{T}^{3}}\left(1+\sum_{r=1}^{T-1}r\alpha(r)\right)$

We additionally impose that $\sum_{r=1}^{\infty}r\alpha(r)<\infty$ because we
use a _decaying function_ for $\alpha(\cdot)$ and $B_{T}^{2}\geq c_{0}T$ with
a positive constant $c_{0}$. We choose

$\alpha(r)\leq C_{1}\frac{1}{r^{2}(\log r)^{p}},$

for $r\geq 2$ and fixed $p>1$, where $0<C_{1}<\infty$ is a constant. Then,

$\displaystyle
T\frac{L^{3}}{B_{T}^{3}}\left(1+\sum_{r=1}^{T-1}r\alpha(r)\right)\leq\underbrace{\frac{L^{3}}{c_{0}^{3/2}}\left(1+\alpha(1)+C_{1}\sum_{r=2}^{\infty}\frac{1}{r(\log
r)^{p}}\right)}_{=K(L,C_{1},c_{0})<\infty}\frac{1}{\sqrt{T}}$ (48)

where $K(L,C_{1},c_{0})$ is a positive and finite constant only depending on
the quantities within the parenthesis. The inequality follows the fact
$B_{T}^{2}\geq c_{0}T$. Thus, we achieve an asymptotic bound of $O(T^{-1/2})$.
This completes our proof. ∎

The weak convergence with $p_{t}$ implies the weak convergence with $q_{t}$,
by a simple application of Lemma 7.2 in Sarkar et al. (2014a).

######  Lemma 7 (Lemma 7.2 in Sarkar et al. (2014a)).

Suppose Lemma 6 holds. Then, under Assumption 1 and assuming $\sigma_{c}>0$,

$\text{Conditioned on }S_{C}\text{, \hskip 10.00002pt
}\sum_{t}q_{t}/\sqrt{T}\overset{d}{\to}\mathcal{N}(0,\sigma_{c}^{2})\text{
\hskip 10.00002pt as }T\rightarrow\infty$

We now prove the weak convergence of our estimator.

###### Proof of Theorem 2.

By the definition of $\widehat{g}_{T}-g$ in Theorem 1, we achieve

$\displaystyle\sqrt{T}(\widetilde{g}_{T}-g)$
$\displaystyle=\sqrt{T}(\widetilde{g}_{T}-\widehat{g}_{T})+\sqrt{T}(\widehat{g}_{T}-{g})$
$\displaystyle=O(\sqrt{T}\beta)+\frac{T^{-1/2}\sum_{t}q_{t}}{\widehat{d}_{T}}+\frac{T^{1/2}B_{T}\mathrm{E}[\widehat{d}_{T}|S_{C}]}{\widehat{d}_{T}}.$

We are able to assess the convergence of each item. Since Lemma 3 shows
$\mathrm{E}[\widehat{d}_{T}|S_{C}]\to R(C)$ and Lemma 2 shows
$\operatorname{Var}(\widehat{d}_{T}|S_{C})\to 0$, we obtain
$\widehat{d}_{T}\overset{p}{\to}R(C)$. By Lemma 7, $T^{-1/2}\sum_{t}q_{t}$
converges to $\mathcal{N}(0,\sigma_{c}^{2})$. Further, Lemma 4 proves
$B_{T}=O(b)$. Combining these results with the Slutsky’s lemma, we get the
following results:

$\displaystyle\sqrt{T}(\widetilde{g}_{T}-g)=O(\sqrt{T}\beta)+\widetilde{W}_{T}+O_{P}(\sqrt{T}b),$

where $\widetilde{W}_{T}$ is a random variable such as

$\displaystyle\widetilde{W}_{T}\overset{d}{\to}\mathcal{N}(0,\sigma_{c}^{2}/R(C)^{2}).$

With the settings $\beta=o(T^{-1/2})$ and $b=o(T^{-1/2})$, we obtain the
statement. ∎

## References

  * Sho (2005) (2005) Chapter 4 particle-particle interactions, in Interfacial Separation of Particles (Eds.) S. Lu, R. J. Pugh and E. Forssberg, Elsevier, vol. 20 of Studies in Interface Science, pp. 172 – 244.
  * Abu-El-Haija et al. (2018) Abu-El-Haija, S., Perozzi, B., Al-Rfou, R. and Alemi, A. A. (2018) Watch your step: Learning node embeddings via graph attention, in Advances in Neural Information Processing Systems, pp. 9180–9190.
  * Adamic and Adar (2001) Adamic, L. A. and Adar, E. (2001) Friends and neighbors on the web, SOCIAL NETWORKS, 25, 211–230.
  * Adamic and Adar (2003) Adamic, L. A. and Adar, E. (2003) Friends and neighbors on the web, Social networks, 25, 211–230.
  * Airoldi et al. (2006) Airoldi, E. M., Blei, D. M., Fienberg, S. E., Xing, E. P. and Jaakkola, T. (2006) Mixed membership stochastic block models for relational data with application to protein-protein interactions, in In Proceedings of the International Biometrics Society Annual Meeting.
  * Barabási and Albert (1999) Barabási, A.-L. and Albert, R. (1999) Emergence of scaling in random networks, science, 286, 509–512.
  * Benson et al. (2018) Benson, A. R., Abebe, R., Schaub, M. T., Jadbabaie, A. and Kleinberg, J. (2018) Simplicial closure and higher-order link prediction, Proceedings of the National Academy of Sciences, 115, E11221–E11230.
  * Bhatia et al. (2019) Bhatia, S., Chatterjee, B., Nathani, D. and Kaul, M. (2019) A persistent homology perspective to the link prediction problem, Springer, vol. 881 of Studies in Computational Intelligence, pp. 27–39.
  * Björner and Kalai (2006) Björner, A. and Kalai, G. (2006) On f-vectors and homology, Annals of the New York Academy of Sciences, 555, 63 – 80.
  * Brin and Page (2012) Brin, S. and Page, L. (2012) Reprint of: The anatomy of a large-scale hypertextual web search engine, Computer networks, 56, 3825–3833.
  * Chen et al. (2010) Chen, L., Goldstein, L. and Shao, Q. (2010) Normal Approximation by Stein’s Method, Springer Verlag.
  * da Xu et al. (2020) da Xu, chuanwei ruan, evren korpeoglu, sushant kumar and kannan achan (2020) Inductive representation learning on temporal graphs, in International Conference on Learning Representations (ICLR).
  * Evans and Harlow (1957) Evans, M. W. and Harlow, F. H. (1957) The particle-in-cell method for hydrodynamic calculations.
  * Greene (2012) Greene, L. H. (2012) Protein structure networks, Briefings in Functional Genomics, 11, 469–478.
  * Grover and Leskovec (2016a) Grover, A. and Leskovec, J. (2016a) node2vec: Scalable feature learning for networks, in Proceedings of the 22nd ACM SIGKDD international conference on Knowledge discovery and data mining, ACM, pp. 855–864.
  * Grover and Leskovec (2016b) Grover, A. and Leskovec, J. (2016b) Node2vec: Scalable feature learning for networks, in Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’16, p. 855–864.
  * Halabi et al. (2009) Halabi, N., Rivoire, O., Leibler, S. and Ranganathan, R. (2009) Protein sectors: evolutionary units of three-dimensional structure, Cell, 138, 774–786.
  * Jeh and Widom (2002) Jeh, G. and Widom, J. (2002) Simrank: a measure of structural-context similarity, in Proceedings of the eighth ACM SIGKDD international conference on Knowledge discovery and data mining, ACM, pp. 538–543.
  * Kipf and Welling (2016) Kipf, T. N. and Welling, M. (2016) Variational graph auto-encoders, NIPS Workshop on Bayesian Deep Learning.
  * Liben-Nowell and Kleinberg (2007a) Liben-Nowell, D. and Kleinberg, J. (2007a) The link-prediction problem for social networks, J. Am. Soc. Inf. Sci. Technol., 58, 1019–1031.
  * Liben-Nowell and Kleinberg (2007b) Liben-Nowell, D. and Kleinberg, J. (2007b) The link-prediction problem for social networks, Journal of the American society for information science and technology, 58, 1019–1031.
  * Lü and Zhou (2011) Lü, L. and Zhou, T. (2011) Link prediction in complex networks: A survey, Physica A: statistical mechanics and its applications, 390, 1150–1170.
  * Mannige (2014) Mannige, R. V. (2014) Dynamic New World: Refining Our View of Protein Structure, Function and Evolution, Proteomes, 2, 128–153.
  * Martinez et al. (2019) Martinez, B., Lobet, M., Duclous, R., d’Humières, E. and Gremillet, L. (2019) High-energy radiation and pair production by Coulomb processes in particle-in-cell simulations, Physics of Plasmas, 26.
  * Mitzenmacher (2004) Mitzenmacher, M. (2004) A brief history of generative models for power law and lognormal distributions, Internet Mathematics, 1, 226–251.
  * Patil et al. (2020) Patil, P., Sharma, G. and Murty, M. N. (2020) Negative sampling for hyperlink prediction in networks, Springer, vol. 12085 of Lecture Notes in Computer Science, pp. 607–619.
  * Perozzi et al. (2014) Perozzi, B., Al-Rfou, R. and Skiena, S. (2014) Deepwalk: Online learning of social representations, in Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining, ACM, pp. 701–710.
  * Pukhov (2015) Pukhov, A. (2015) Particle-in-cell codes for plasma-based particle acceleration.
  * Rossi et al. (2020) Rossi, E., Chamberlain, B., Frasca, F., Eynard, D., Monti, F. and Bronstein, M. (2020) Temporal graph networks for deep learning on dynamic graphs, in ICML 2020 Workshop on Graph Representation Learning.
  * Salton and McGill (1986) Salton, G. and McGill, M. J. (1986) Introduction to Modern Information Retrieval, McGraw-Hill, Inc., USA.
  * Sarkar et al. (2014a) Sarkar, P., Chakrabarti, D. and Jordan, M. (2014a) Nonparametric link prediction in large scale dynamic networks, Electron. J. Statist., 8, 2022–2065.
  * Sarkar et al. (2014b) Sarkar, P., Chakrabarti, D., Jordan, M. et al. (2014b) Nonparametric link prediction in large scale dynamic networks, Electronic Journal of Statistics, 8, 2022–2065.
  * Shalaby et al. (2017) Shalaby, M., Broderick, A. E., Chang, P., Pfrommer, C., Lamberts, A. and Puchwein, E. (2017) SHARP: A spatially higher-order, relativistic particle-in-cell code, The Astrophysical Journal, 841, 52.
  * Sharma et al. (2020) Sharma, C., Chauhan, J. and Kaul, M. (2020) Learning representations using spectral-biased random walks on graphs.
  * Shervashidze et al. (2011) Shervashidze, N., Schweitzer, P., van Leeuwen, E. J., Mehlhorn, K. and Borgwardt, K. M. (2011) Weisfeiler-lehman graph kernels., J. Mach. Learn. Res., 12, 2539–2561.
  * Stein et al. (1986) Stein, C., (Organization), J. and of Mathematical Statistics, I. (1986) Approximate Computation of Expectations, IMS Lecture Notes.
  * Sunklodas (2007) Sunklodas, J. (2007) On normal approximation for strongly mixing random variables, Acta Applicandae Mathematicae, 97, 251–260.
  * Van der Vaart (2000) Van der Vaart, A. W. (2000) Asymptotic statistics, vol. 3, Cambridge university press.
  * Villani (2008) Villani, C. (2008) Optimal transport: old and new, vol. 338, Springer Science & Business Media.
  * Wegscheider (1911) Wegscheider, R. (1911) Über simultane gleichgewichte und die beziehungen zwischen thermodynamik und reactionskinetik homogener systeme, Monatshefte für Chemie und verwandte Teile anderer Wissenschaften, 32, 849–906.
  * Xu et al. (2013) Xu, Y., Rockmore, D. and Kleinbaum, A. M. (2013) Hyperlink prediction in hypernetworks using latent social features, in Discovery Science (Eds.) J. Fürnkranz, E. Hüllermeier and T. Higuchi, Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 324–339.
  * Yoon et al. (2020) Yoon, S., Song, H., Shin, K. and Yi, Y. (2020) How much and when do we need higher-order informationin hypergraphs? a case study on hyperedge prediction, Proceedings of The Web Conference 2020.
  * Zhang and Chen (2018a) Zhang, M. and Chen, Y. (2018a) Link prediction based on graph neural networks, in Advances in Neural Information Processing Systems, pp. 5165–5175.
  * Zhang and Chen (2018b) Zhang, M. and Chen, Y. (2018b) Link prediction based on graph neural networks, arXiv preprint arXiv:1802.09691.
  * Zhang et al. (2018) Zhang, M., Cui, Z., Jiang, S. and Chen, Y. (2018) Beyond link prediction: Predicting hyperlinks in adjacency space.
  * Zhou et al. (2009a) Zhou, T., Lü, L. and Zhang, Y.-C. (2009a) Predicting missing links via local information, The European Physical Journal B, 71, 623–630.
  * Zhou et al. (2009b) Zhou, T., Lü, L. and Zhang, Y.-C. (2009b) Predicting missing links via local information, The European Physical Journal B: Condensed Matter and Complex Systems, 71, 623–630.