Year: 2020

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 higherorder prediction methods. As a theoretical achievement, we prove the consistency and asymptotic normality in terms of Wasserstein distance of our estimator using Stein’s method.

heading: 1 INTRODUCTION
text: Numerous types of networks like social (Liben-Nowell & 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′ > 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 1 (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 ht at time t, the inference task is to predict the future arrival of a new hyperedge ht′ , which covers a larger set of vertices than ht and contains all the vertices in ht. Figure 1 (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 two major 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 & Adar, 2003; Brin & Page, 2012; Jeh & Widom, 2002; Zhou et al., 2009a; Barabási & 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 1 (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). This problem is remedied by the use of simplices because they naturally model these substructures as a collection of subsets (i.e., faces) of the simplex.
We provide real-world examples of where our simplicial complex based approach can play a significant role. (i) Organic Chemistry: It is quite common to have the same set of elements interacting with each other in different configurations, which result in very different functioning compounds (Ma et al., 2011). Specifically, R-thalidomide and S-thalidomide are two different configurations of thalidomide, where the R-form was meant to help sedate pregnant women, while the S-form unfortunately resulted in birth defects. This is a famous example in stereo chemistry to show the consequences of mistaking two extremely close configurations (differing by a single bond) as being the same. Structure prediction to avoid such phenomenon in drug synthesis allows chemists to achieve a much higher yield and avoid wastage of expensive resources. (ii) Gene expression networks: Gene networks have nodes that represent genes and edges connect genes which have similar expression patterns (Zhang & Horvath, 2005). Subgraphs called modules are tightly connected genes in such a gene expression network. Genomics research provides evidence that higher-order gene expression relationships (like second and third-order) and their measurements can have very important implications for cancer prognosis. When making structural predictions in these aforementioned examples, our simplicial complex based approach provides much more fine-grained control over competing methods by capturing subtler differences in configurations.
To combat these challenges, our approach views the evolving graph1 as a time process under the framework of nonparametric time series prediction, which models the evolution of higher-order structures (as simplices) and their local neighborhoods (spatial dimension) over a moving time window (temporal dimension). Our inference problem is then modeled as predicting the evolution to a higher-dimensional simplex at time t′ > t, given a simplex at time t. It is important to note that this task is more general and greatly diverges from the task proposed by Benson et al. (2018). Our
1We 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.
task requires just a single simplex σ in order to predict a higher-dimensional simplex τ , whose face / subset is σ, whereas Benson et al. (2018) requires the presence of all constituent σ faces in order to predict τ . For example, in Benson et al. (2018) (also shown in Figure 1 (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 1 (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].
To this effect, we succinctly capture the features characterizing the local neighborhood of a simplex as a combination of a face-vector (Björner & Kalai, 2006), which is a well-established vector signature in combinatorial topology literature and a novel scoring function which infers the affinity of sub-simplices based on the strength of their past interactions. Based on these features, we design a kernel estimator to infer future evolution to higher-dimensional simplices and prove both the consistency and asymptotic normality of our estimator.
Our contributions: (a) We propose a kernel estimator that predicts higher-order interactions in an evolving network. (b) We prove the consistency and asymptotic normality of our kernel estimator. (c) We evaluate our method on real-world dynamic networks by proposing higher-order link prediction baselines and observe significant gains in prediction accuracy in comparison to the baselines.

heading: 1.1 RELATED STUDIES
text: 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 & Adar, 2003), PageRank (Brin & Page, 2012), SimRank (Jeh & Widom, 2002), resource allocation (Zhou et al., 2009a), preferential attachment (Barabási & Albert, 1999), persistence homology based ranking (Bhatia et al., 2019), and similarity-based methods (Liben-Nowell & Kleinberg, 2007b; Lü & Zhou, 2011). (ii) Random walk based methods consist of DeepWalk (Perozzi et al., 2014) , Node2Vec (Grover & Leskovec, 2016b) 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 & Welling, 2016), WYS (Abu-El-Haija et al., 2018), and SEAL (Zhang & Chen, 2018b).
Higher-order link prediction: 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). 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 graph simplicial complex is closed under taking subsets, which enables us to better encode more information for improved inference.

heading: 2 PRELIMINARY: GRAPH SIMPLICIAL COMPLEX (GSC)
text: 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 σ is an element of A, then so is every nonempty subset of σ. The element σ 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) G on G is an ASC consisting of subsets of V . In particular, G is a collection of subgraphs of G. With graphs, we denote a ddimensional simplex (or d-simplex) of a GSC by σ(d) = [v0, v1, . . . , vd]. Each non-empty subset of σ(d) is called a face of σ(d).
We define several notions related to GSCs, that are useful for describing the evolution of graphs. Definition 2 (Filtered GSC). For I ⊂ N, a filtered GSC indexed over I is a family (Gt)t∈I of GSCs such that for every t ≤ t′ in I , Gt ⊂ Gt′ holds.
Obviously, Gt0 ⊂ Gt1 ⊂ . . .Gtn is a discrete filtration induced by the arrival times of the simplices: Gti \ Gti−1 = σti . 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 Gt,p for the last p discrete time steps is defined as Gt,p := (Gt′)tt′=t−p = (Gt−p, . . . ,Gt), where Gt′ ⊃ Gt′−1.
We define a notion for dealing with the neighborhood around a given simplex σ(d) ∈ G. We introduce a set of all simplices of dimension d′ or less from G, i.e., G(d ′) − := {σ(d) ∈ G | d ≤ d′}. G (0) − is a vertex set and G(1)− is the set of edges and vertices. We also write i ∼ j whenever vertices i and j are adjacent in G(1)− , and write i ∼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 G(1)− . Then, we define a ball around σ(d). Definition 3 (k-ball centered at vertex and simplex). At time t, we define a k-ball centered at vertex i by Bt,k(i) := {j : i ∼k j and i, j ∈ G(0)t− }. We also define a k-ball centered at a simplex σ(d) as Bt,k(σ (d)) := ⋃ i:Vert(σ(d))Bt,k(i), where Vert(σ (d)) denotes all vertexes in σ(d).
Now, we define a sub-complex G′t(σ(d)) ⊆ Gt as the GSC that contains all the simplices in Gt spanned by the vertices in the k-ball Bt,k(σ(d)).

heading: 3 PREDICTING THE ARRIVAL OF HIGHER-ORDER SIMPLICES
text: We consider the prediction of a simplex’s arrival in the setting described below. Consider a filtered GSC Gt,p. At time t, given a d-dimensional simplex σ(d) = [v0, · · · , vd], we predict the formation of a (d + 1)-simplex τ (d+1) = [v0, · · · , vd, ṽ] with a new vertex ṽ ∈ Bt,k(σ(d)). Here, we restrict ṽ to be k-reachable from σ(d). In order to find out which simplices and vertices are most likely to appear in τ (d+1) at time t+ 1, we need to design features for σ(d) and ṽ.

heading: 3.1 FEATURE DESIGN FOR SIMPLEX
text: 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 invariant2 of the GSC. Definition 4 (face-vector). A combinatorial statistic of G is the face-vector (or f -vector) as f(G) = (f−1, f0, . . . , fd−1), where fk = fk(G) records the number of k-dimensional faces σ(k) ∈ G.
We then define the feature of a simplex σ(d) at time t, denoted by Nt(σ(d)) = f(G′t(σ(d))). In words, the feature is compactly represented as the face-vector of sub-complex G′t(σ(d)). Our facevector 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 ṽ using its proximity to σ(d). To this end, we begin by describing affinity between two vertices. Given two vertices
2A topological invariant is a property that is preserved by homeomorphisms.
v, v′ ∈ G(0)− , we denote by s(v, v′) the weighted sum of all past co-occurrences of vertices v and v′ in σ(d), where the weight is d from σ(d). We then devise a scoring function h(·, ·) that assigns an integral score to the possible introduction of a vertex ṽ to a d-simplex σ(d) = [v0, · · · , vd] as
ht(σ (d), ṽ) = d∑ i=0 s(vi, ṽ). (1)
Intuitively, the scoring function describes higher co-occurrence between σ(d) and ṽ at time t, indicates a higher likelihood of forming a (d + 1)-simplex τ (d+1) = [σ(d), 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 σ(d) at time t, and a possible introduction of a new vertex ṽ ∈ Bt,k(σ(d)), we assign a feature vector
Ft(σ (d), ṽ) = (Nt(σ (d)), ht(σ (d), ṽ)). (2)
We denote the set of all such possible (σ(d), ṽ) pairs with their corresponding feature vectors equal to F as Pt(σ(d), F ). Furthermore, among the pairs in Pt(σ(d), F ), we denote by P τt (σ
(d), F ) those set of pairs with feature vectors equal to F that actually form τ (d+1) = [σ(d), ṽ] at time t, i.e., σ(d) appears as a face in τ (d+1) at time t. Note the distinction that not all d-simplices counted in Pt(σ(d), F ) end up being promoted to higher (d + 1)-simplices in the next time step. Table 1 provides a list of notations.

heading: 3.2 PREDICTION MODEL AND KERNEL ESTIMATOR
text: For the prediction, we define an indicator variable that displays the appearance of a new simplex. Given a d-simplex σ(d) = [v0, · · · , vd] ∈ Gt, the arrival at time t+ 1 of a (d+ 1)-simplex τ (d+1) = [v0, · · · , vd, ṽ] with a new vertex ṽ ∈ Bt,k(σ(d)) is captured by the following variable
Yt+1(τ (d+1)) := { 1 if σ(d) is a face of τ (d+1) , 0 otherwise.
(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:
Yt+1(τ (d+1)) | Gt,p ∼ Bernoulli(g(Ft(σ(d), ṽ))) (4)
where 0 ≤ g(·) ≤ 1 is a function of the feature vector Ft(σ(d), ṽ). In words, the indicator variable Yt+1(τ (d+1)) is Bernoulli distributed with success probability given by function g(Ft(σ(d), ṽ)), conditioned on having seen the last p states of GSC Gt. This model describes that the appearance probabilities for two simplices σi and σ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 Gt(d) be a set of d-dimensional simplices at time t from Gt,p. Also, for brevity, let F represent a feature Ft(σ(d), ṽm) with some t, σ(d) and ṽm subject to ṽm ∈ Bt,k(σ(d)i ). Let ‖F − F ′‖1 denote the L1-distance between two feature vectors, and also define a L1-ball Γ(F, δ) := {F ′ : ‖F − F ′‖1 ≤ δ}. We define our kernel function K(·, ·) as follows. With two features F and F ′, we define it as
K (F, F ′) := I {F = F ′}+ βI {‖F − F ′‖1 ≤ δ}
1 + β|Γ(F, δ)| , (5)
where β > 0 is the bandwidth parameter, and 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 L1-ball of radius δ 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 σ(d) and ṽ and set F = FT (σ(d), ṽ), our estimator of g(F ) in (Equation 4) is written as follows:
g̃T (F ) :=
∑T t′=T−p ∑ σ (d) j ∈Gt′ (d) ∑ ṽn∈Bt′,k(σ (d) j ) K ( F, Ft′(σ (d) j , ṽn) ) · Yt′+1([σ(d)j , ṽn])∑T
t′=T−p ∑ σ (d) j ∈Gt′ (d) ∑ ṽn∈Bt′,k(σ (d) j ) K ( F, Ft′(σ (d) j , ṽn) ) . (6) Time-complexity of our estimator: As both the face-vectors and counts Pt+1 / P (τ) t+1 (updated in a data cube) are computed simultaneously, they incur the same time overhead. At time t, for a dsimplex σ(d) it takes O(|Vk| + |Ek|) time to compute a k-ball around σ(d), where Vk and Ek are the set of vertices and edges in the k-hop subgraph of a vertex. We must check this k-ball of σ(d) against p|G(d)t− | 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(|Vk|) time. Recall that fd(Gt) denoted the total number of d-simplices in Gt. Vk and Ek are the set of vertices and edges in the k-hop subgraph of a vertex. We must check this k-ball of σ(d) against p|G(d)t− | 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(Tpfd(Gt)(|Vk|+ |Ek|))|G(d)t− |). Storage-complexity of our estimator: For a single simplex, our estimator requires storing: (i) a pair of integer counts, namely (Pt(·, ·), P τt (·, ·)), in a datacube, which costs O(1) and (ii) a d + 1- dimensional face vector which takes storage O(d+ 1). As the total number of simplices is denoted by fd(Gt), we arrive at a total storage cost of O(dfd(Gt)).
4 THEORETICAL PROPERTY OF THE ESTIMATOR
𝑔𝑔"!
𝑇 → ∞
𝑔𝑔"! Consistency Asymptotic Normality We show that our estimator has theoretical validity: (i) consistency and (ii) asymptotic normality. (i) The consistency guarantees g̃T achieves zero error as T increases by converging to g. (ii) The asymptotic normality implies that the error 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).

heading: 4.1 CONSISTENCY
text: 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 S. We denote the time of entering class C by TC and the event as E(TC). Let SC denote the event St ∈ C, where St is the state of the Markov chain at time t. Then, E(TC) ∩ SC is the event that the chain enters class C at time TC and remains in that communication class indefinitely.
With the event SC , we provide the bias-variance decomposition of 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 ĥT (F ) = (T − p)−1 ∑T−1 t=p ∑|Gt(d)| j=1 |Gt(d)| −1|P τt+1(σ (d) j , F )| and d̂T (F ) = (T −
p)−1 ∑T−1 t=p ∑|Gt(d)| j=1 |Gt(d)|
−1|Pt+1(σ(d)j , F )|. For the sake of brevity, we fix F and denote g̃T (F ) by g̃T . Similarly, g, ĥT and d̂T are used. Also, we define a term BT (F,C) = E[ĥT | SC ]/E[d̂T | SC ]− g. Then, by Proposition 2 in the supplementary material, we decompose g̃T − g as
g̃T − g = [ĥT − gd̂T ]− E[ĥT − gd̂T | SC ]
d̂T︸ ︷︷ ︸ =:VT (variance)
+ BT (F,C)E[d̂T | SC ]
d̂T︸ ︷︷ ︸ =:BT (bias)
+o(1). (7)
To bound the variance term VT , we make an assumption that our Markov chain Xt exhibits a α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 α-mixing: Definition 5 (α-mixing). A stochastic process Xt is α-mixing, if a coefficient α(r), defined as
α(r) = sup |t1−t2|≥r
{|Pr(A ∩B)− Pr(A)Pr(B)| : A ∈ Σ(X−t1), B ∈ Σ(X + t2)},
satisfies α(r) → 0 as r → ∞. Here, Σ(X−t1) and Σ(X + t2) are the sigma-algebras of past and future events of the stochastic process up to and including t1.
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 BT , 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 κ : R → R in the Schwartz space (i.e. it is infinitely differentiable and converging to zero faster than any polynomial as x→ ±∞) with a parameter b > 0 such that |g(F )− g(F ′)| = O(κ(−‖F − F ′‖1/b)), as b→ 0, ∀F, F ′.
We utilize all the assumptions, then prove the consistency of g̃T . Theorem 1 (Consistency). Suppose that the GSC filtration process is α-mixing, β = o(1), and Assumption 1 holds. Then, for any F and conditional on SC , our estimator g̃T (F ) is well-defined with probability tending to 1, and |g̃T (F )− g(F )| p→ 0 holds as T →∞.

heading: 4.2 ASYMPTOTIC NORMALITY
text: 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 Zn by a Gaussian random variable, where Zn is a sum of n mean-centered random variables {Ai}ni=0, where Ai corresponds to a random variable which depends on the i-th d-simplex σ(d)i in our GSC, which is dependent on other d-simplices
whose neighborhoods largely overlap with that of σ(d)i . Here, let dw 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 Zn. We provide the following theoretical result. Its formal statement is Theorem 3, which is deferred to the supplementary material due to its complexity.
Proposition 1 (Gaussian approximation for dependent variables; Simple version of Theorem 3). Suppose Zn = ∑n i=0Ai, where {Ai}ni=0 is generated by zero-mean random variables
X0, X1, · · · , Xn satisfying the α-mixing condition, such as Ai = Xi/Bn with Bn = E[ ∑n i=0X 2 i ]. Also, suppose that Pr{|Xi| ≤ L} = 1 holds for i = 1, · · · , n with some constant L > 0. Then, with an existing finite constant C > 0, we have dw(Zn, N) ≤ C (∑n i=1E|Ai|3 + TL3
B3n
∑n−1 i=1 nα(n) ) .
This result extends Sunklodas (2007) in the sense of Markov chains on GSCs that satisfy the α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.
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 VT and BT from (7). Then, we regard VT as a sum of dependent random variables and apply the result developed in Proposition 1. Let σ2c be a limit of variance of the numerator in T−1/2VT as T → ∞. We recall that SC denotes the event St ∈ C, where St is the state of the Markov chain at time t. Theorem 2 (Asymptotic Normality). Suppose that Assumption 1 holds, the GSC filtration process is α-mixing, and σc > 0. If β = o(T−1/2) and b = o(T−1/2), then, for any F and conditioned on SC , the following holds: √ T (g̃T (F )− g(F )) d→ N (0, σ2c/R(C)2), as T →∞.
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.

heading: 5 REAL-WORLD DATA EXPERIMENTS
text: 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′ > t, which contains σ(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.8K, x = 5K), EU (n = 998, m = 29.3K, x = 8K), Contact (n = 327, m = 5.8K, x = 10K), and NDC (National Drug Code) (n = 1.1K, m = 6.2K, x = 12K).
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 β, 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 & Adar, 2001), Jaccard Coefficient (JC) (Salton & McGill, 1986), and Preferential attachment (PA) (Mitzenmacher, 2004), (ii) deep-learning based (Node2vec (NV) (Grover & Leskovec, 2016a), and SEAL (SL) (Zhang & Chen, 2018a)), 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.

heading: 5.1 RESULTS AND DISCUSSION
text: 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 β 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.
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 3, 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 (σ(d), v) pair formation depending on the dimension of the simplex in which the pair co-occur in the past.

heading: 6 CONCLUSION
text: 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.

heading: A EXAMPLE OF SIMPLEX AND RELATED NOTIONS
text: Example 1. We begin by computing the k-balls centered at 1-simplex [9, 10] in Gt and Gt−1, respectively.
The k-ball at time t for k = 1 (i.e., 1-hop vertices only) centered at [9, 10] is:
Bt,1([9, 10]) = Bt,1([9]) ∪ Bt,1([10]). This is the union of k-balls at underlying vertices 9 and 10 according to Definition 3.
Bt,1([9, 10]) = Bt,1([9]) ∪Bt,1([10]) = {9, 13, 8, 5, 6, 10} ∪ {10, 14, 13, 9, 6, 7, 11, 15} = {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:
Bt−1,1([9, 10]) = Bt−1,1([9]) ∪Bt−1,1([10]) = {9, 13, 8, 5, 6, 10} ∪ {10, 14, 13, 9, 6, 11, 15} = {9, 13, 8, 5, 6, 10, 14, 11, 15}
Notice that there is only a difference of vertex 7 missing from set Bt−1,1([9, 10]) as compared to set Bt,1([9, 10]) at time t.
Now, we calculate the subcomplex spanned by Bt,1([9, 10]) Then3,
G′t([9, 10]) ={[9], [13], [8], [5], [6], [10], [14], [11], [15], [7], (8) [9, 10], [9, 6], [9, 5], [9, 8], [9, 13], [10, 6], [10, 7], [10, 11], [10, 15], [10, 14], [10, 13],
[9, 6, 10], [10, 14, 15]}
Now, we compute the compressed f -vector notation of G′t([9, 10]) to get
f(G′t([9, 10])) = (1, 10, 11, 2)
Finally, the neighborhood Nt([9, 10]) = (1, 10, 11, 2).

heading: B FURTHER DETAILS OF THE EXPERIMENT
text: B.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.8K, x = 5K) and EU (n = 998, m = 29.3K, x = 8K) 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.8K, x = 10K) 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.1K, m = 6.2K, x = 12K) 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.
B.2 COMPARED METHODS
Adamic-Adar Adamic & Adar (2001) and the Jaccard Coefficient Salton & 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
3all simplices are placed on a line each in increasing order of their dimension
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 & Leskovec (2016a) and SEAL Zhang & Chen (2018a) 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 sayA andB at time t, we predict whether authorsA andB will co-author with a third author C (ternary co-authorship relationship) on a single paper at time t+ 1 in the future.

heading: C PROOF FOR CONSISTENCY
text: For preparation, we rewrite the estimator ĝT . Plugging in the definition of our kernel (Equation 5) into the equation of our estimator (Equation 6) along with the definitions of Pt(·, ·) and P τt (·, ·) to replace the indicator variables with actual counts, we obtain the following simplification of Equation 6. Then, it is reformulated as
g̃T (F ) =
∑T t′=T−p ∑ σ (d) j ∈Gt′ (d) ( |P τt′+1(σ (d) j , F )|+ β ∑ s∈Γ(F,δ)|P τt′+1(σ (d) j , s)| ) ∑T t′=T−p ∑ σ (d) j ∈Gt′ (d) ( |Pt′+1(σ(d)j , F )|+ β ∑ s∈Γ(F,δ)|Pt′+1(σ (d) j , s)|
) . (9) When we set β = 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′ + 1. This fraction is summed across various d-simplices σ(d)j by varying j and also across various discrete time steps by varying t′. Setting β > 0 allows our estimator to smooth over close by features.
It turns out simpler to study a proxy estimator ĝT , which omits the smoothing of the original one. We will show that it is asymptotically equivalent to g̃. Let |Gt(d)| denote the total number of d-simplices in GSC Gt at time t. Then, for a feature F , we define the proxy estimator as
ĝT (F ) := ĥT (F )
d̂T (F ) , (10)
where the terms are defined as
ĥT (F ) = 1
T − p T−1∑ t=p |Gt(d)|∑ j=1 |P τt+1(σ (d) j , F )| |Gt(d)| , and d̂T (F ) = 1 T − p T−1∑ t=p |Gt(d)|∑ j=1 |Pt+1(σ(d)j , F )| |Gt(d)| .
Observe that ∑|Gt(d)| i=1 ranges over the total number of d-simplices in Gt4 and note that |Gt(d)| changes with time step t. Recall that the terms P τt+1(·, ·) and Pt+1(·, ·) count all actual and possible formation of [σ(d)i , ṽm] ((d + 1)-simplex) that have the same feature vector F . Lemma 1 in the supplementary material proves |g̃T (F ) − ĝT (F )| → 0 as β → 0. First of all, we show the validity of the proxy estimator ĝT .
Lemma 1 (Approximation by proxy). We obtain
|g̃T (F )− ĝT (F )| = O(β),∀F. 4In practice, the total number of d-simplices is much less than the maximum possible cliques with n vertices.
Proof of Lemma 1. Recall that Γ(F, δ) denotes the set of features at L1-distance at most δ from F . We denote by |Γ(F, δ)| the cardinality of this set.
g̃T (F ) = ĥT (F ) + C1
d̂T (F ) + C2 where C1 = β ∑ j,t ∑ s∈Γ(F,δ)|P τt+1(σ (d) j , s)| and C2 = β ∑ j,t ∑ s∈Γ(F,δ)|Pt+1(σ (d) j , s)|. Due to the finiteness of features in P τt+1(·, ·) and Pt+1(·, ·), we have that C1 = C2 = O(β). Both C1 and C2 are non-negative integers. So,
|g̃T (F )− ĝT (F )| = ∣∣∣∣∣ ĥT (F ) + C1d̂T (F ) + C2 − ĥT (F )d̂T (F ) ∣∣∣∣∣ = O(β)
In the last step, the second fraction is a positive constant and can thus be ignored from our asymptotic analysis because both ĥT and d̂T are bounded.
Next, we prove the convergence of the proxy estimator 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
ĝT − g = VT + BT .
Furthermore, with the event SC , there exist a stochastic terms qt for t such as
VT = (T − p)−1
∑T−1 t=p qt
d̂T .
Proof of Proposition 2. With the definition of BT (F,C), we have
ĝT − g = ĥT
d̂T − g (11)
= ĥT − gd̂T
d̂T
= [ĥT − gd̂T ]− E[ĥT − gd̂T | SC ] + E[ĥT − gd̂T | SC ]
d̂T
= [ĥT − gd̂T ]− E[ĥT − gd̂T | SC ] d̂T + E[ĥT | SC ]− gE[d̂T | SC ] d̂T
= [ĥT − gd̂T ]− E[ĥT − gd̂T | SC ] d̂T + BT (F,C)E[d̂T | SC ]
d̂T = VT + BT . (12)
We are interested in the asymptotic behavior of the Markov chain at time T → ∞. Let F denote FT (σ (d) i , ṽm). Recall that the terms |P τT+1(·, ·)| and |PT+1(·, ·)| count all actual and possible formation of [σ(d)i , 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 SC for any communication class C.
For a given time step t, we define
ĥT (t) := 1
|Gt(d)| |Gt(d)|∑ j=1 |P τt+1(σ (d) j , F )| (13)
d̂T (t) := 1
|Gt(d)| |Gt(d)|∑ j=1 |Pt+1(σ(d)j , F )|
Note that ĥT = 1T−p ∑T−1 t=p ĥT (t) and d̂T = 1 T−p ∑T−1 t=p d̂T (t).
Let us set
qt := [ĥT (t)− gd̂T (t)]− E[ĥT (t)− gd̂T (t) | SC ] (14)
Note that qt is the numerator of the stochastic term in Equation 7 and a bounded deterministic function of SC at a given time step t.
For the stochastic term 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 α-mixing, then, as T →∞, we obtain
Var(ĥT (F )|SC)→ 0, and Var(d̂T (F )|SC)→ 0,
for any F .
Proof of Lemma 2. We show that variance divided by T converges to a non-negative constant. Let UT := ∑ t qt/ √ T , where qt (as shown in Equation 14) is a bounded deterministic function of the state of Xt at time t. As demonstrated in Sarkar et. al. Sarkar et al. (2014a), we too break our weighted sum UT across three time intervals: (i) [1, TC − 1], (ii) [TC , TC + M − 1], and (iii) [TC +M,T ], where M is a constant. Now, from Sarkar et al. (2014a), we simply apply Lemma 5.7 to get that E[Var(UT | E(TC), SC) | SC ] → σc, for some σc ≥ 0 and from Lemma 5.8, we have that Var(E[UT | E(TC), SC ] | SC) = o(1). Now, since the law of total variance provides
Var(UT | SC) = E[Var(UT | E(TC), SC) | SC ] + Var(E[UT | E(TC), SC ] | SC)
we use the previous results from Lemmas 5.7 and 5.8 in Sarkar et al. (2014a) to get
Var(UT | SC)→ σc as T → 0, for some constant σc ≥ 0
Plugging in the definition of qt into UT and calculating Var(UT | SC), it follows trivially that Var(ĥT | SC)→ 0 and Var(d̂T | SC)→ 0 as T →∞. 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 α-mixing, then there exist a function R(C) with a deterministic function of class C denote, such as
lim T→∞ E[d̂T (F ) | E(TC), SC ] = R(C), and lim T→∞ E[d̂T (F ) | SC ] = R(C).
Proof of Lemma 3. We know by definition that
E[d̂T (F ) | E(TC), SC ] = 1
T − p T−1∑ t=p |Gt(d)|∑ j=1 E [ |Pt+1(σ(d)j , F )| |Gt(d)| ∣∣∣∣ E(TC), SC ]
(15)
This is an average of termsE [ |Pt+1(σ(d)j ,F )| |Gt(d)| ∣∣∣∣ E(TC), SC] spanning across d-simplices with indices j ∈ {1, · · · , |Gt(d)|} and discrete time steps t ∈ {p, · · · , T − 1}. For ease of notation, let
Xj := |Pt+1(σ(d)j , F )| |Gt(d)|
Xj denotes the total number of possible (d+ 1)-simplices with a d-face as σ (d) j divided by the total number of d-simplices in Gt.
In the R.H.S. of Equation 15, the term inside the summation is simplified as E[Xj | E(TC), SC ] = ∑ x xPr[Xj = x | E(TC), SC ]
We know that both Pt+1(σ (d) j , F ) and Gt(d) are fully determined given the current state St of the Markov chain. Let IS(Y ) denote an indicator variable of whether “Y is in state S” or not.
We have, Pr[Xj = x | E(TC), SC ] = ∑ S IS(Xj = x)Pr[St = S | E(TC), SC ]
As a result, the R.H.S. of Equation 15 becomes
1
T ∑ t ∑ S ( ∑ j,x xIS(Xj = x))Pr[St = S | E(TC), SC ] (16)
Let λ(S) = ∑ j,x xIS(Xj = x) as this term is fully determined by state S. Then, Equation 16 can be rewritten as ∑ S λ(S) ∑ t Pr[St = S | E(TC), SC ] T (17)
Due to stationarity, the average ∑ t Pr[St = S | E(TC), SC ]/T will converge to a constant function of state S, denoted by R(S). Given that λ(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 limT→∞E[d̂(F ) | SC ] = RC . 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 α-mixing, then, conditional on SC , we obtain VT p→ 0 as T →∞.
Proof of Proposition 3. By Proposition 2, the term VT is written as (T−p)−1
∑ t qt
d̂T .
For the denominator, Lemma 3 shows that E[d̂T | SC ] → R(C), where R(C) is a positive deterministic function of class C. Also, Lemma 2 states that V ar(d̂T | SC)→ 0 holds as T →∞. Thus, d̂T
p→ R(C) > 0 holds. Also, VT is asymptotically well defined for class C. For the nominator, Lemma 2 also shows that limT→∞ V ar( ∑ t qt/T | SC) = 0 as T → ∞ and
E[qt | SC ] = 0, therefore we have limT→∞ 1/T ∑ t qt
qm→ 0 conditioned on SC . By the continuous mapping theorem, we obtain the statement.
Next, we discuss the bias term BT . To this aim, we rewrite the term BT (F,C) as follows:
BT (F,C) = (E[ĥT (F ) | SC ]− g(F )E[d̂T (F ) | SC ])
E[d̂T (F ) | SC ] .
Lemma 4. Given that Assumption 1 holds and that when T →∞, the bandwidth parameter b→ 0. Then, we have that BT (F,C) = O(b) = o(1) as T →∞.
Proof of Lemma 4. For t ∈ [p, T − 1], j ∈ [1, |Gt(d)|] and a fixed feature vector F , the numerator of BT (F,C) can be expressed as an average of terms of the form
At := E
[ |P τt+1(σ (d) j , F )|
|Gt(d)|
∣∣∣∣ SC ] − E [ |Pt+1(σ(d)j , F )| |Gt(d)| ∣∣∣∣ SC ] g(F ) (18)
The first term in Equation 18 can be rewritten using the tower property as
E [ E [ |P τt+1(σ (d) j , F )|
|Gt(d)|
∣∣∣∣ E(TC), SC ] ∣∣∣∣ SC ]
When we condition on E(TC), it makes |P τt+1(σ (d) j ,F )|
|Gt(d)| conditionally independent of SC , if t > TC . Also, for t ≥ TC , we have
E
[ |P τt+1(σ (d) j , F )|
|Gt(d)|
∣∣∣∣ E(TC), SC ] = |Pt+1(σ(d)j , F )| |Gt(d)| g(Ft(σ (d) j , ṽn)) (19)
where ṽn ∈ Bt−1,k(σ(d)j ).
Given the result in Equation 19 and the fact that the term |P τt+1(σ (d) j ,F )|
|Gt(d)| is bounded results in
E
[ |P τt+1(σ (d) j , F )|
|Gt(d)|
∣∣∣∣ E(TC), SC ]
(20)
≤ |Pt+1(σ(d)j , F )| |Gt(d)| g(Ft(σ (d) j , ṽn))I[TC ≤ t] + cI[TC > t]
≤ |Pt+1(σ(d)j , F )| |Gt(d)| g(Ft(σ (d) j , ṽn)) + cI[TC > t],
where c > 0 an existing constant.
Now, the numerator of BT (F,C) can be upper bounded as∣∣∣∣∣∑ t At/T ∣∣∣∣∣ ≤∑ t 1 T ∣∣∣∣∣E [ |Pt+1(σ(d)j , F )| |Gt(d)| (g(Ft(σ (d) j , ṽn))− g(F )) ∣∣∣∣ SC ]∣∣∣∣∣ (21)
+ c′ ∑ t Pr[TC > t]/T.
The second term in Equation 21 vanishes as T →∞ because it is of order O(E[TC ]/T ). Thus, the numerator of BT (F,C) is an average of terms of the form
E [ |Pt+1(σ(d)j , F )| |Gt(d)| (g(Ft(σ (d) j , ṽn))− g(F )) ∣∣∣∣ SC ] . (22)
Our feature vector counts and simplex neighborhoods are finite because |Gt(d)| is bounded. The expectation in Equation 22 is just a summation of finite terms. We set F ′ = Ft(σ (d) j , ṽn), and make use of our smoothness assumption 1, so that
|g(F ′)− g(F )| = O(κ(−‖F − F ′‖1/b)). We also use Lemma 3 to say that the denominator of our bias term converges to a constant R(C). So,
BT (F,C) = O(κ(−‖F − F ′‖1/b)) = O(b). The last equality follows the property of κ in the Schwartz space. Now, since BT (F,C) = O(b) and b→ 0 as T →∞, then we have that BT (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 SC , BT p→ 0 holds as T →∞.
Proof of Proposition 4. Proposition 2 shows that BT (F ) = BT (F,C)/d̂T (F ).
By Lemma 2 and 3, we obtain d̂T → R(C) > 0 as T → ∞, as similarly shown in the proof of Proposition 3. Further, Lemma 4 states that BT (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 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 g̃T , we additionally combine the result of Lemma 1, then obtain the statement.

heading: D PROOF FOR ASYMPTOTIC NORMALITY
text: D.1 INTRODUCTION TO WASSERSTEIN DISTANCE AND APPROXIMATION TECHNIQUE
We denote by BL(R) the space of such bounded functions h that are 1-Lipschitz. More formally,
‖h‖∞ = sup x∈R |h(x)| <∞ and Lip(h) = 1, where Lip(h) = sup x 6=y |h(x)− h(y)| |x− y|
So, h ∈ BL(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
dw(W,Z) = sup h∈BL(R)
|Eh(W )− Eh(Z)| (23)
D.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
Ef ′(W ) = E[Wf(W )], (24)
for all absolutely continuous functions f : R → R with E|f ′(Z)| < ∞. Conversely, if Equation 24 holds for all bounded, continuous, and piecewise continuously differentiable functions f with E|f ′(Z)| <∞, 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 : R→ R. Here, Stein’s lemma (Lemma 5) shows that W d= Z, if
Ef ′(W )− 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
f ′(W )−Wf(W ) = h(W )− Eh(Z) (26)
The f that satisfy Equation 26 with h ∈ BL(R) must satisfy the following conditions for all y, z ∈ R
‖f‖ ≤ 2 , ‖f ′‖ ≤ 2 , ‖f ′′‖ ≤ √ 2/π (27)
|f ′(y + z)− f ′(y)| ≤ D|z|
where h0(y) = h(y)−Eh(Z), c1 = supx≥0 ξ(x), c2 = supx≥0 x(1−xξ(x)), and ξ(x) = (1−Φ)/φ (where Φ(x) is the distribution function and φ(x) = Φ′(x)). Then, D = (c1 + c2)‖h0‖∞ + 2 is
a constant. Additionally, we have a bound on the covariance, given the dependent r.v.’s are also bounded.
If Pr{|X| ≤ C1} = Pr{|Y | ≤ C2} = 1, then (28) |Cov(X,Y )| ≤ 4C1C2α(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 α 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).
D.2 GAUSSIAN APPROXIMATION FOR DEPENDENT VARIABLES WITH GSC
In our model, we assume the r.v. Ai to represent a d-simplex σ (d) i in a GSC G. In order to have a notion of α-mixing in our setting, we must first define a distance between two d-simplices σi and σj . We drop the (d) superscript for brevity and ease of notation. We define this distance as the Hausdorff distance between simplices as
dH(σi, σj) = max { sup v∈σi inf v′∈σj dg(v, v ′), sup v′∈σj inf v∈σi dg(v, v ′) } (29)
where dg(v, v′) counts the number of edges in the geodesic connecting vertices v and v′ in G(1)− . In Stein’s method, the sum of dependent r.v.’s is studied by breaking the sum Zn into two sets based on the r in mixing coefficient α(r). In our setting, given a fixed d-simplex σi, we study two partial sums pertaining to: 1) all d-simplices that are at most r-apart from Ai and 2) the remaining partial sum after removing the variables pertaining to 1) from Zn.
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 α-mixing in Markov chains based on GSCs. For a sequence of r.v.’s X1, X2, · · · satisfying the α-mixing condition, we write
Zn = n∑ i=1 Ai, Ai = Xi Bn , B2n = E( n∑ i=0 Xi) 2
We assume Bn > 0.
T (m) i denotes the contribution of d-simplices that are further than m away from σi and x(σi, r)
denotes the partial sum of r.v.’s representing simplices that are exactly r away from Ai. Therefore,∑m r=0 x(σi, r) gets us all those d-simplices that are greater than or equal to m away from σi. We are interested in the contribution of simplices r away from σi as we vary r from 0 to m. With Proposition 5, we proceed to derive an upper bound on dw(Zn, N) (i.e., the Wasserstein distance between Zn and N ). Proposition 5. Let S(σi, r) denote the set of d-simplices whose Hausdorff distance equals r. More formally, S(σi, r) = {σj : dH(σi, σj) = r} Additionally, let X̂ denote a mean-centered version of r.v. X . Then,
x(σi, r) = ∑
p∈S(σi,r)
Ap (x(σi, 0) = Ai)
T (m) i = Zn − m∑ r=0 x(σi, r), m = 0, 1, · · · , (T (−1)i = Zn)
Suppose that EZn = 0, EZ2n = 1, and EA 2 i < ∞ for all i = 1, · · · , n. Let be a r.v. uniformly distributed in [0, 1] and independent of other r.v.’s. Let f : R→ R be a differentiable function such that supx∈R |f ′(x)| <∞. Then we have
Ef ′(Zn)− EZnf(Zn) = E1 + · · ·+ E7
where E1 = − n∑ i=1 ∑ r≥1 EAix(σi, r) [ f ′(T (r) i + x(σi, r))− f ′(T (r) i ) ] E2 = −
n∑ i=1 EA2i [ f ′(T (0) i + Ai)− f ′(T (0) i ) ] E3 = −
n∑ i=1 ∑ r≥1 2r∑ q=r+1 E ̂Aix(σi, r)δ(q)i , E4 = − n∑ i=1 ∑ r≥1 ∑ q≥2r+1 E ̂Aix(σi, r)δ(q)i
E5 = n∑ i=1 ∑ r≥1 EAix(σi, r) r∑ q=0 Eδ (q) i , E6 = − n∑ i=1 ∑ q≥1 EÂ2i δ (q) i , E7 = n∑ i=0 EA2iEδ (q) i ,
and
δ (q) i = f ′(T (q−1) i )− f ′(T (q) i ).
Theorem 3. Consider a sequence of r.v.’sX1, X2, · · · that satisfy α-mixing condition (Definition 5). Let EXi = 0, Pr{|Xi| ≤ L} = 1, for i = 1, · · · , n, for some constant L > 0. Then, for every h ∈ BL(R),
dw(Zn, N) ≤ C(S,D) ( n∑ i=0 E|Ai|3 + nL3 B3n n−1∑ r=1 rα(r) ) where C(S,D) is a finite constant which depends on D and maxi |S(σi, r)|.
Proof of Theorem 3. Given a d-simplex σi and its corresponding r.v. Ai, recall that S(σi, r) denotes the set of d-simplices that are at Hausdorff distance r away from σi. We additionally define Sm to denote the maximum cardinality of S(σi, r) for all i and a fixed r, i.e.,
Sm = max i |S(σi, r)|
In order to upper bound the Wasserstein distance between Zn and N , we estimate the difference Eh(Zn) − Eh(N) using Proposition 5. It was shown in Proposition 5 that this difference is a sum of terms E1, · · · , E7. We will proceed by individually bounding each term. Bounding E1:
|E1| = ∣∣∣∣∣∣∣∣ n∑ i=1 ∑ r≥1 EAix(σi, r) [ f ′(T (r) i + x(σi, r))− f ′(T (r) i ) ] ︸ ︷︷ ︸
(i) ∣∣∣∣∣∣∣∣ (30) We can upper bound term (i) in Equation 30 using Equation 27 by
D| ||x(σi, r)| ≤ D|x(σi, r)| (since | | ≤ 1) Then, ∣∣∣x(σi, r) [f ′(T (r)i + x(σi, r))− f ′(T (r)i )]∣∣∣ ≤ D(x(σi, r))2 (31) Now, we upper bound x(σi, r) as
x(σi, r) ≤ Sm ( L
Bn
) (32)
since each normalized r.v. in Zn is upper bounded by L/Bn. The L.H.S. of Equation 31 is upper bounded by DS2m
( L2
B2n
) .
We know that ∣∣∣∣∣Cov ( Ai︸︷︷︸, x(σi, r) [f ′(T (r)i + x(σi, r))− f ′(T (r)i )]︸ ︷︷ ︸ )∣∣∣∣∣ (33)
= E [ Aix(σi, r) [ f ′(T (r) i + x(σi, r))− f ′(T (r) i ) ]] ︸ ︷︷ ︸
(ii)
+ EAi︸︷︷︸ =0
Ex(σi, r) [ f ′(T (r) i + x(σi, r))− f ′(T (r) i ) ]
Notice that term (ii) is nothing but the summand in Equation 30. We apply the covariance bounds (Equation 28), to obtain
(ii) ≤ 4 ( L
Bn
)( DS2m L2
B2n
) α(r) ≤ 4DS2m L3
B3n α(r) (34)
Then, we have that
|E1| ≤ n∑ i=1 ∑ r≥1 ( 4DS2m L3 B3n α(r) ) ≤ 4DS2m nL3 B3n n−1∑ r=1 α(r). (35)
Here, the summation ∑n−1 r=1 appears, because there are only n variables that should measure the dependence between each other.
Bounding E2:
|E2| = ∣∣∣∣∣∣∣∣ n∑ i=1 EA2i [ f ′(T (0) i + Ai)− f ′(T (0) i ) ] ︸ ︷︷ ︸
(iii) ∣∣∣∣∣∣∣∣ (36) Using Equation 27, we have that term (iii) ≤ D|Ai|. Then, |E2| can simply be bounded as
|E2| ≤ D n∑ i=1 E|Ai|3 (37)
Bounding E3:
|E3| = ∣∣∣∣∣∣ n∑ i=1 ∑ r≥1 2r∑ q=r+1 E ̂Aix(σi, r)δ(q)i ∣∣∣∣∣∣ (38) ≤
n∑ i=1 ∑ r≥1 2r∑ q=r+1 ∣∣∣EAix(σi, r)δ(q)i ∣∣∣︸ ︷︷ ︸ (iv) + n∑ i=1 ∑ r≥1 2r∑ q=r+1 |EAix(σi, r)|E ∣∣∣δ(q)i ∣∣∣︸ ︷︷ ︸ (v)
Let us focus on bounding terms (iv) and (v), separately.
For term (iv): We can further split it as∣∣∣∣∣∣∣E Ai︸︷︷︸a x(σi, r)δ (q) i︸ ︷︷ ︸
b ∣∣∣∣∣∣∣ (39) We have previously worked out the bounds for terms Ai and x(σi, r). Therefore, we now focus our attention on bounding δ(q)i .
δ (q) i = f ′(T (q−1) i )− f ′(T (q) i ) (40)
= f ′(T (q) i + x(σi, r))− f ′(T (q) i ) ≤ D|x(σi, r)| (using Equation 27)
≤ D ( Sm L
Bn
)
Therefore, term b (in term (iv)) is upper bounded by DS2m L2
B2n .
Applying the covariance bound (Equation 28), we have that (iv) ≤ 4 ( L
Bn
)( DS2m L2
B2n
) α(r) ≤ 4DS2m L3
B3n α(r) (41)
For term (v): We have calculated some bounds previously, so we can again split (v) as∣∣∣∣E Ai︸︷︷︸x(σi, r)︸ ︷︷ ︸ ∣∣∣∣︸ ︷︷ ︸
c
E ∣∣∣δ(q)i ∣∣∣︸ ︷︷ ︸ d ≤ [ 4 ( L Bn )( Sm L Bn ) α(r) ] ︸ ︷︷ ︸
using Eqn 28 on c
[ D ( Sm L
Bn )] ︸ ︷︷ ︸
using Eqn 40 on d
(42)
≤ 4DS2m L3
B3n α(r)
Now, combining the inequalities for terms (iv) and (v), we have that
|E3| ≤ n∑ i=1 ∑ r≥1 2r∑ q=r+1 8DS2m L3 B3n α(r) ≤ 8DS2m nL3 B3n n−1∑ r=1 rα(r) (43)
Bounding E4: We have that
|E4| = ∣∣∣∣∣∣ n∑ i=1 ∑ r≥1 ∑ q≥2r+1 E ̂Aix(σi, r)δ(q)i ∣∣∣∣∣∣ (44) ≤ 4DS2m nL3
B3n ∑ r≥1 ∑ q≥2r+1 α(q − r)
≤ 4DS2m nL3
B3n n−1∑ r=1 rα(r)
because ∑ r≥1 ∑ q≥2r+1 α(q − r) < ∑n−1 r=1 rα(r).
Now, using our previous bounds, the remaining terms, i.e., |E5|, |E6|, and |E7| are bounded as follows:
|E5| ≤ 4DS2m nL3
B3n n−1∑ r=1 rα(r) (45)
|E6| ≤ 4DS2m nL3
B3n n−1∑ r=1 α(r)
|E7| ≤ D n∑ i=0 E|Ai|3
Finally, we combine these bounds to arrive at our final upper bound as
|Ef ′(Zn)− EZnf(Zn)| ≤ C(Sm, D) ( n∑ i=0 E|Ai|3 + nL3 B3n n−1∑ r=1 rα(r) ) (46)
where C(Sm, D) is a constant term depending on constants Km and D. This completes our proof.
D.3 DEFERRED PROOF
In this section, we establish our estimator’s result on the asymptotic normality.
Let n = |Gt(d)| denote the total number of d-simplices in GSC Gt at time t. Recall from Equation 13 that
ĥT (t) := 1
n n∑ j=1 |P τt+1(σ (d) j , F )|
d̂T (t) := 1
n n∑ j=1 |Pt+1(σ(d)j , F )|
Also, from Equation 14, we have
qt := [ĥT (t)− gd̂T (t)]− E[ĥT (t)− gd̂T (t) | SC ]
Additionally, let us define a variable pt for convenience as follows
pt := [ĥT (t)− gd̂T (t)]− E[ĥT (t)− gd̂T (t) | E(TC), SC ]
Note that pt is conditioned on both E(TC) and SC , while qt is just conditioned on SC . Keeping these expressions in mind, we begin by showing the weak convergence with pt. Lemma 6. Under Assumption 1, given pt for any finite TC and σc > 0∑
t≥TC+M
pt/ √ T d→ N (0, σ2c ) conditioned on E(TC) ∩ SC as T →∞
Proof of Lemma 6. Given a normalized r.v. WT which is a sum of weakly dependent r.v.’s as
WT :=
∑ t≥TC+M pt√
Var (∑ t≥TC+M pt | E(TC), SC )
where pt 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
WT d→ N (0, 1) conditioned on E(TC) ∩ SC
Recall that by Lemma 5.7 in Sarkar et al. (2014a), we have
Var  ∑ t≥TC+M pt | E(TC), SC  /T → σ2c Now, we show that conditioned on event E(TC) ∩ SC , our bound in Theorem 3 has a convergence rate of O(T−1/2).
Note that our pt corresponds toAi and T to n in Theorem 3. As pt is a function of St, it involves p+1 GSCs (Gt−p+1, · · · ,Gt+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(σi, r) only 2 states that are at distance r apart from σi, i.e., Ai − r and Ai + r. Therefore, Sm = O(1). Given that Pr{|pt| ≤ L} = 1, for t = 1, · · · , T and for some constant bound L > 0, we have
T∑ t=1 E|pt|3 + T L3 B3T T−1∑ r=1 rα(r) (ignoring some constant terms) (47)
≤ T L 3
B3T
( 1 +
T−1∑ r=1 rα(r) ) We additionally impose that ∑∞ r=1 rα(r) < ∞ because we use a decaying function for α(·) and B2T ≥ c0T with a positive constant c0. We choose
α(r) ≤ C1 1
r2(log r)p ,
for r ≥ 2 and fixed p > 1, where 0 < C1 <∞ is a constant. Then,
T L3
B3T
( 1 +
T−1∑ r=1 rα(r)
) ≤ L 3
c 3/2 0
( 1 + α(1) + C1
∞∑ r=2
1
r(log r)p ) ︸ ︷︷ ︸
=K(L,C1,c0)<∞
1√ T
(48)
where K(L,C1, c0) is a positive and finite constant only depending on the quantities within the parenthesis. The inequality follows the fact B2T ≥ c0T . Thus, we achieve an asymptotic bound of O(T−1/2). This completes our proof.
The weak convergence with pt implies the weak convergence with qt, 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 σc > 0,
Conditioned on SC , ∑ t qt/ √ T d→ N (0, σ2c ) as T →∞
We now prove the weak convergence of our estimator.
Proof of Theorem 2. By the definition of ĝT − g in Theorem 1, we achieve √ T (g̃T − g) = √ T (g̃T − ĝT ) + √ T (ĝT − g)
= O( √ Tβ) +
T−1/2 ∑ t qt
d̂T + T 1/2BTE[d̂T |SC ] d̂T .
We are able to assess the convergence of each item. Since Lemma 3 shows E[d̂T |SC ]→ R(C) and Lemma 2 shows Var(d̂T |SC) → 0, we obtain d̂T p→ R(C). By Lemma 7, T−1/2 ∑ t qt converges to N (0, σ2c ). Further, Lemma 4 proves BT = O(b). Combining these results with the Slutsky’s lemma, we get the following results:
√ T (g̃T − g) = O( √ Tβ) + W̃T +OP ( √ Tb),
where W̃T is a random variable such as
W̃T d→ N (0, σ2c/R(C)2).
With the settings β = o(T−1/2) and b = o(T−1/2), we obtain the statement.

Review summaries:
This paper presents an estimator that predict higher-order structure in time-varying graphs. The authors present an kernel-based estimator, prove that it is consistent when the indicator variable for whether a particular (d+1)-dimensional simplex is Bernoulli distributed with a function g. The authors prove that their estimator is asymptotically normal.  The authors also present some experiments on real-world data

This paper provide a method for high-order structure prediction problem. Specifically, the paper first defines a high-order structure on graphs named graph simplicial complex (GSC). Then the paper introduces a feature generation method used for the high-order structures. The features are also used in the proposed method for high-order structure prediction. The proposed method is based on a nonparametric kernel which carries the feature similarities of high-order structures. With this kernel the method uses a Bernoulli distribution for the prediction of the existence of the high-order structure in unseen times.