Intensity Profile Projection: A Framework for Continuous-Time Representation Learning for Dynamic Networks
Authors: Alexander Modell
Ian Gallagher
Emma Ceccherini
Nick Whiteley
Patrick Rubin-Delanchy
Venue: NeurIPS 2023 poster (2023)
Keywords: dynamic networks
representation learning
spectral methods
TL;DR: We present a new algorithmic framework, Intensity Profile Projection, for learning continuously-evolving representations of the nodes of a dynamic network.

Abstract: We present a new representation learning framework, Intensity Profile Projection, for continuous-time dynamic network data. Given triples $(i,j,t)$, each representing a time-stamped ($t$) interaction between two entities ($i,j$), our procedure returns a continuous-time trajectory for each node, representing its behaviour over time. The framework consists of three stages: estimating pairwise intensity functions, e.g. via kernel smoothing; learning a projection which minimises a notion of intensity reconstruction error; and constructing evolving node representations via the learned projection. The trajectories satisfy two properties, known as structural and temporal coherence, which we see as fundamental for reliable inference. Moreoever, we develop estimation theory providing tight control on the error of any estimated trajectory, indicating that the representations could even be used in quite noise-sensitive follow-on analyses. The theory also elucidates the role of smoothing as a bias-variance trade-off, and shows how we can reduce the level of smoothing as the signal-to-noise ratio increases on account of the algorithm `borrowing strength' across the network.

Intensity Profile Projection: A Framework for Continuous-Time Representation Learning for Dynamic Networks
Alexander Modell
Ian Gallagher Emma Ceccherini Nick Whiteley?
Patrick Rubin-Delanchy"
'Imeeial College London, U.K. University of Bristol, U.K.
a. modell @imperial aaliaall........ ac. uk, {ian. galagher,, emma. ceccherini, nick. whiteley, matrick.rubinnna.. aawaanycyyyyyrroe ac
Abstract
We present new representation learning framework, Intensity Profile Projection, for continuous-time dynamic network data. Given triples (i, J, (), each representing a time-stamped (I) interaction between two entities (), our procedure returns a continuous-tin trajectory for each node, representing its behaviour over time. The framework consists of three stages: estimating pairwise intensity functions, e.g. via kernel smoothing; learning projection which minimises notion of intensity reconstruction error; and constructing evolving node representations via the Iearned projection. The trajectories satisfy two properties, known structural and temporal coherence, which we see as fundamental for reliable inference. Moreoever, we develop estimation theory providing tight control on the error of any estimated trajectory, indicating that the representations could even be used in quite noise- sensitive follow-on analyses. The theory also elucidates the role of smoothing bias-variance trade-off, and shows how we can reduce the level of smoothing as the signal-to-noise ratio increases account of the algorithm 'borrowing strength' across the network.
Introduction
Making sense of patterns of connections occurring over time is common theme of modern data analysis and is often approached in one two ways. On the one hand, we may see dynamic network data as graph, in which connections between the same entities over time are somehow treated as one, e.g. through weighting. This view evokes methodological ideas such as community detection (12, topological data analysis (3), manifold learning 4 ( On the other, may the data as a set of point processes [6]. each modelling the event times of connections between two entities. This view evokes temporal notions such trend, changepoints and periodicity.
this paper, we develop a representation learning framework for continuous- time dynamic network data in which ideas from both the graph and temporal domains can be combined. The framework obtains a continuously evolving trajectory for each node which represents continuously evolving behaviour in the network. These trajectories tion tasks such as spatio-temporal clustering, behavioural analytics, detecting network trends and periodicities, and forecasting. Our framework provides two basic guarantees which could reasonably be viewed as minimum requirements for these tasks:
be used for data exploration, inference and predic-
structural coherence: when two nodes exhibit similar behaviour at a point time, their representations at that time are close;
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

2. temporal coherence: when node exhibits similar behaviours at two points in time, its representations those times are close.
Despite this, almost all existing embedding algorithms fail to satisfy both properties, and there is a perception that some trade-off between the two is necessary 7 [. To our knowledge, unfolded spectral embedding
[ 10 is the only exception, but the method is limited to discrete time.
At its core, our procedure is a spectral method and can be implemented at scale using a sparse singular value decomposition. We give preliminary motivation for the algorithm as minimising notion of reconstruction error, before developing more rigorous estimation theory: a non-asymptotic. uniform error bound under minimal assumptions the data generating process, which draws on recent developments in entrywise eigenvector estimation for random matrices [10]44]. This tight control of the behaviour of estimated trajectories means that even quite noise-sensitive follow-on analyses, such as topological data analysis, could plausibly be shown to be consistent.
With appropriate development, the ideas of this paper could be brought to bear on several problems. The first is building a coherent narrative about the nodes based on an embedding. In contact-tracing data for interactions between children at school (Section [5), we see trajectories mixing during lunch hours and then returning to their earlier positions, which represent physical classrooms, Both types of coherence are necessary for this description to be possible. We cannot see this, for example, just by looking at the pattern of clustering over time. Given richer and larger datasets it could be possible to automate this "story-telling" process. Among many possible applications, such technology could transform how we use contact-tracing data in the future. Second, the jump from discrete continuous time has substantial implications for mathematical modelling, since makes notions such as continuity and smoothness possible. This work could pave the way to describing structural dynamics in formal mathematical and quantitative terms, for example, what mean by 'communities splitting" (15), including the exact time at which this occurs, or by "network polarisation" (16), and how could be measured, e.g. as a type of derivative. Finally, in criminal investigations, dynamic networks are often analysed to locate an entity, such as a victim of human- trafficking, whose pseudonym or other identifier changing [17). Similarly, in dynamic networks of corporate contracting, it is useful to detect when one company is acting like another in the past, e.g. when a shell company takes over the illegal operations of a sanctioned company [18). Since two nodes acting the same way, at different epochs, are still embedded to the same position, our framework could lead to novel matching and tracking technologies.
Related work. While modelling and performing inference on continuous-time dynamic networks has a well-established literature [6] [1). the majority of existing methods for representation learning obtain a single, static representation of each node (22-30). and we only aware of two existing methods which learn continuously evolving node representations from the data we consider [3]] 32]
Representation learning for discrete-time dynamic networks has a more established literature, and algorithms broadly fall under community detection methods [3), fitting latent position models 37 spectral methods p 10
42 and word-embedding-- bassddddgg---ddd methods 43 44.
For specific
choice of intensity estimator (the histogram), our method can be viewed as weighted graph analogue of unfolded spectral embedding p 10], a spectral method for discrete-time and multilayer networks, but those papers consider different data and models. Given how limited the options are for handling continuous time, our method comparison we also include some discrete-time methods which could reasonably be used as alternatives.
2 Intensity Profile Projection
Data. We consider dynamic network data, denoted G. representing instantancous undirected interac- tions between nodes over time, which we define formally as G = (V. E) on a time domain T = (0, T), containing a vertex set V = [n]] and a set of triples& & = {(Te. je. te (e)eee each corresponding to an undirected interaction event, where ie < je E V. te E T. We let Egg == {t (i.j.t) E E denote the interaction events between nodes and j.
Model. We assume the interaction events Eig are driven by an independent inhomogeneous Poisson process with intensity (t) Informally:
(t)dt =P (interaction between nodes and j in (t, f + dt|}

---

Algorithm Approximate Intensity Profile Projection
Input: Continuous time dynamic graph G, dimension
1: Construct intensity estimates (²5(-) of A-5(.) from Eyg for all i <.
2: Compute the top d left singular vectors Ui...., Udd of
((t))) A(t2)
A(tB)
t equally spaced points (0, 7) 3: Define the trajectory node as
where
x= = (
where Üa := (T,,.... ... Qd). Output: Node trajectories {,,))....,.(((,...
We represent these intensities in symmetric time-varying matrix A(-) T R Procedure. Intensity Profile Projection can be summarised follows.
1. Intensity estimation. Construct intensity estimates )-)(-) of )sy(-) from Egy for all i < j.
2. Subspace learning. Compute the top d cigenvectors Ü. = (ăI...., (2) of
Σ 3 7 // (2)d),,
where (t) has symmetric entries Rıg (t), and rows denoted (,(1) called intensity profiles.
3. Projection. For query node at time t, project the intensity profile (t( onto the subspace spanned by "1...., Ga to obtain x2(() (.A((()
While we develop more principled statistical justifications for the procedure in future sections, is inspired by simple reconstruction argument. For an arbitrary d'dimensional subspace spanned by the orthonormal columns of matrix Va E Rxd let
| VvV - R.(t))l
(t: Va)
denote the reconstruction error of node time and define the integrated residual sum of squares as
(V(va) = T EE (l.Vaad
Lemma 1. Among all d-dimensional subspaces of R", the column span of U minimises the integrated residual sum of squares criterion Ria
Lemma D may be viewed as dynamic analogue the classical Eckart- Young theorem on low-rank matrix approximation [45). A proof is given in Section of the appendix.
2.1 Intensity estimation
The choice of intensity estimator is left fully open, but our theory makes two important recommenda- tions. First, there are computational gains to be made using sparse estimators for subspace learning. Second, the procedure borrows strength across the network, and can give precise representations even when the individual intensity estimates are noisy (e.g. inconsistent). In our experiments, we focus on standard non oo-aammmtriic estimators such the histogram kernel smoothers, and choose kernels with finite support to induce sparse estimates.
2.2 Subspace learning
The subspace learning step of our procedure involves the computation of an integral, and computing the cigendecomposition of the resulting dense matrix Σ, both of which may be infcasible for large

---

networks. If sparse intensity estimator is employed in step I of the procedure and we approximate the integral ( using numerical quadrature scheme, then step 2 can rephrased as a single sparse, truncated singular value decomposition, which can be computed quickly for very large networks using a efficient solver [46 47
Consider the numerical approximation
A (tb)
where t <.. < tg are equally spaced points on (0, T]. The top d eigenvectors of the right-hand hand- side of (2) are then equal" the top left singular vectors of the matrix
A(ti) Att2) A(tB))
(3)
the row concatenation of (4)...... ((B)) This procedure is presented in Algorithm and we discuss its computational complexity in Section Df the appendix.
2.3 Projection
The inductive nature of the Intensity Profile Projection allows to obtain representations (t) demand, for example, the full trajectory for particular node, or the representations of the entire graph at point in time. It is possible obtain representations for intensity profiles outside the training sample, corresponding to new nodes or times outside the training domain, allowing online inference. In practice, the discussion of this computational and statistical trade-off for future work (see, for example, [88 in the context of static networks).
will need retrain occasionally, i.e. return step 2, although we leave
Estimation theory
In this section, we develop estimation theory showing the sense in which {t(t) "good" estimator of Xt) U (t) where Ua (u1,..., Ud) E Raxd is the matrix containing the top-d
orthonormal eigenvectors
J (t)ddt.
E = 7
In this section, we assume, without loss of generality, that T= = quantities which appear our main theorem. Firstly, assume that each Xxy()) Lipschitz with constant L, and is upper bounded by Amax- Secondly, we define the (reduced) condition number and the eigengap.
(0, 1]. We now introduce some
=
and
== od dd+1s
σ
respectively, where of 2 o eigenvalues of Σ. Finally, introduce the subspace coherence parameter
U 22
which is small when, informally, information about single entry of is "spread out" across the matrix 44
Rather than attempt develop theoretical framework encompassing all intensity estimators, choose arguably the most rudimentary, the histogram, and we expect more powerful estimators will only improve matters. This choice of estimator also attractive because allows us pinpoint the crucial practical considerations at play.
Up to signs, rotations in the eigenspaces in the case of repeated eigenvalues. and assuming gap between the dth and (d 1)th eigenvalues.

---

Notation. We say an event E occurs with overwhelming probability if P(E) 2 -I-n" for any constant c > We use S denote that a v Cb where > is universal constant which, when qualified with the prior probabilistic statement, may depend on the constant c. Additionally, we write a b if a S b and a 2
We now state the assumptions we require for our theorem. Our first assumption is that the intensities are bounded.
Assumption (Bounded intensities). The intensities are upper bounded by a constant which doesn't 't depend on the other quantities in the problem; i.e. Amax X 1.
Our second assumption is on the population integrated residuals. It ensures that the intensity profiles A (4),..... (1) do not deviate "too much" from a common low-dimensional subspace. Assumption 2 (Small population residuals). The population residuals satisfy
d (t())....,r((() & V μ log"/2
for all t E [0, 1].
Our third assumption is a technical condition on the eigengap which, broadly speaking, ensures that there is "enough signal".
Assumption (Enough signal). The eigengap satisfies å log(ó/ 0/
/ /TAmax)) Z Kn.maxx
Our final assumption is on the bin size, and ensures that the bins are not chosen "too small".
Assumption 4 (Large enough bins). The number of bins satisfies M S n1max/ log n.
These assumptions weaker than those typically required in the literature (e.g. stochastic block models and random dot product graphs (2 (50). To emphasise this point, consider the following stronger alternative assumptions which imply Assumptions I and 3
Assumption 1a. The intensities (1) are of comparable order, i.e. (1)(() = p for some p S 1 and (0, 1].
all i.j [m], E
Assumption 2a. The matrix E has rank d x 1; is incoherent, i.e. / 1; and its non-zero eigenvectors are of comparable order, i.e. σi x 0a > da+1 0.
is immediate that Assumption Ta implies Assumption D and under Assumption 2a the population residuals are all exactly zero, K and np, which implies Assumptions 2 and 3
Assumption requires that the expected number events involving each node in each bin is at least of the order log" n. This is analogous to the log n degree growth required for perfect clustering under the binary stochastic block model. Since the latter is an information-theores bound [5] and the additional logarithmic powers in our work stem from the sub-exponential tails of the Poisson distribution, we do not think this assumption can be weakened.
We now state our main theorem, which under Assumptions 4 provides non-asymptotic bound on the error between the learned representations and their population counterparts, which holds uniformly over the whole node-set and the time domain.
Theorem 1. Suppose that dig(t) are histogram estimates with M equally-spaced bins and that Assumptions 74 hold. Then with overwhelming probability, there exists an orthogonal matrix W such that
T/3/2 LAmax + Pv /M Ammxxd 3/2 n.
-x(0())x((())
max sup xE|m] [E(0.1]
M8
The proof of Theorem 1 is given in Section of the appendix. As corollary to Theorem I
state a simplified version of this result in which we replace Assumptions 1 and 3 with the stronger Assumptions Ta and 2a Since the Lipschitz constant I scales with the order of the intensities, and we define the quantity Lo satisfying L = pLD which is invariant to the rescaling of intensities.
Corollary 1. Suppose that (t)(1) are histogram estimates with M equally-spaced bins and that Assumptions Ta 2a and hold.. Then with overwhelming probability, there exists an orthogonal matrix W such that
n /22ppL
max sup sEjn 1E (.X,(t) x,(t) S
((x())))x((t)))
+ / Mp log n.
(5)
M
5

---

bins
bins
bins
blas
6
10
12
2
4
6
10
12
2
4
6
10
12
variance
6
12
10 12
time
time
X,(t)
8((t)
8,1t)
Figure bias-variance trade-off. We simulate a network with common intensities (t) 0.7 x {2 + cos(t)} for all i.j. and apply Intensity Profile Projection with histogram intensity estimator with 5, 20, and 200 bins. In the 'bias' plots, the gray lines shows an estimand (t).,, while the blue lines shows its histogram approximation. The discrepancy between the gray line and the blue line corresponds the bias of the Intensity Profile Projection estimator. the blues lines
the 'variance' plots, in the 'bias' plots and the orange line show the estimate obtains using Intensity Profile Projection into one dimension. The discrepancy between the blue line and the orange line corresponds the variance of the Intensity Profile Projection estimator.
3.1
A birinnariance trade-off
The first term in the bound corresponds to the bias between X;(() and X,((). where Xi(() is histogram approximation to Xt(t) (modulo orthogonal transformation, see Section F of the appendix).
The second term corresponds to the variance of the estimate.
Theorem D gives some theoretical guidance on how select the number of bins in the histogram estimator. For simplicity, we consider the setting of Corollary [ bound (5) is optimised by choosing
Ignoring logarithmic terms in n, the
M (mpL3) 1/3
Figure illlustrate this bias- trade-off with example. We simulate a dynamic network with 100 nodes with common intensities Ai (t) = 0.7 x (2 + (ts( for all i,j, the time domain 4T.
The top row shows the population representation X²(t) of single node (gray) and its histogram approximation Xx(() (blue) for a variety of bin sizes. The more bins that are chosen, the smaller the bias and the more (t) resembles (t). The bottom rows shows the histogram approximation (t), and the estimate (t) (orange) obtained using Intensity Profile Projection. The fewer bins that are chosen, the smaller the variance and the more that Xt(t) resembles {t)t).
4 Structural and temporal coherence
For many practical inference tasks, is desirable for representation learning procedure to possess the following two properties:
Structural coherence. two nodes exhibit statistically indistinguishable behaviour at a given time, then their representations at that time are close. That if A,(t) At),,, then
Xi(1)= xt):);
Temporal coherence. If node exhibits statistically indistinguishable behaviour two distinct points in time, then its representations at both these times are close. That is. if (s(s) At), then (s) x (t)..

---

(a) IPP (kernel smoothing)
(b) IPP (histogram) USE
(c) Aligned
(d) Average
0.2 0.4
0.8
0.2 0.4 0.6
0.8
0.2 0.4 0.6 time
0.2
0.8
time
time
time
(e) CLPM
(f) Omnibus
(g) PisCES
(h) MultiNeSS
0.2
0.2
0.2
0.2
time
time
time
time
Community
Communiry?
Figure 2: One-dimensional PCA visualisation of the two-dimensional node representations, obtained using a collection of methods, for a network simulated from a bifurcating block model. Colours correspond to the community membership of the node.
has been observed in a recent survey of 7 that almost all existing dynamic network embedding procedures possess only one of these properties, but not both.
In the following lemma, formally define X,(s) = (t) to mean that (s)s) (t),,, referring to equality "up to statistical noise" in the sense of Theorem [
Lemma 2. Intensity Profile Projection is both structurally and temporally coherent.
This follows from the simple observation that X!( is fixed function of A,(() for all i and To the best of our knowledge, Intensity Profile Projection only existing continuous-time procedure which satisfies these desiderata.
4.1 Simulated example: bifurcating block model
To illustrate these properties, we simulate a two-community dynamic stochastic block model (i.e. where (t) is block structured) in which the nntracoommunity intensities and inter-community intensities initially distinct, they then gradually merge, remain indistinguishable for some time, and finally diverge. We refer to this model as a bifurcating block model and provide full details the simulation in Section of the appendix.
We apply Intensity Profile Projection to the simulated network. using both histogram intensity estimator, and a kernel smoother, to produce two-dimensional representations. For visualisation, we reduce the dimension from two to one using a dynamic adaptation of principal component analysis (see Section of the appendix), and the resulting representations are shown in Figures 2[a) and (b).
In both cases, the estimated trajectories mirror the underlying dynamics of the network: the two communities are in well separated to begin with, gradually merge, remain relatively constant before returning the positions in which they started.
We now illustrate the potential pitfalls of some more naive approaches for embedding dynamic networks. We find that most existing methodology can be viewed as some combination of the two techniques:
Obtain a sequence of static snapshots the network, embed each of the networks snapshots separately and subsequently align the embedding from window + with the embedding from window t.
Alignment 44
Averaging 25 Obtain static summary of the network by averaging it over time, and to embed this to obtain constant node representations.
7

---

Alignment is structurally coherent, however can fail to be temporally coherent. Averaging is tem- porally coherent, but can fail be structurally coherent. To illustrate this point, we apply both approaches, using adjacency spectral embedding into two dimensions, orthogonal Procrustes align- ment and linear interpolation, to a network simulated from the bifureating block model. Figures c) and (d) show visualisations of the trajectories obtained from each approach.
4.2 Method comparison
In this section, we demonstrate how our procedure compares some existing methods on the simulated data described above. Due to the limited number of continuous- time methods, we include number of discrete-time methods (Omnibus, PisCSE and MultiNeSS) which we give an input discrete sequence of snapshots ( (MM of our simulated continuous-time networks. We compare the following methods:
IPP (kernel smoothing). Algorithm | applied with intensities estimated using kernel smoothing.
IPP (histogram) USE [T]]. Algorithm applied with intensities estimated using D
histogram estimator: Equivalent to weighted extension of the Unfolded Spectral Embedding algorithm of [0].
CLPM 31. Fits continuous latent position model log (t) (z,((t (t)l
with a penalty on large velocities the latent space.
Omnibus 40). Approximately factorises the matrix A with blocks A [k, = (Aa)
+
A()),, using spectral decomposition.
PisCES 41. Minimises the objective function
M [L()) L*))ll +
M-1
[L' - ((00) (*
+ 1)
Km
for L'(1), .... (L), LMM,,, where E [o. 1] and L(k) are the Laplacian normalisations of A (k). Then, approximately factorises each L (1), (I) (M) using spectral decompositions.
MultiNeSS [52]. Fits a latent position model (ky(k) (() Φ), where Q(:0, Φ) is a parametric distribution.
We use an embedding dimension of d = 2 for all methods, and for visualisation we reduce this to one using PCA. Additional details such as hyperparameter selection, where applicable, are given in the Section of the appendix.
The CLPM and Omnibus methods produce representations which are temporally coherent, how- ever both fail to capture the complete merging of the communities, shown by Figures ( and (f), and are therefore not structurally coherent. The PisCES and MultiNeSS methods produce repre- sentations which are structurally coherent, however both are unstable when the communities are indistinguisaabl, shown by Figures and (h), and are therefore not temporally coherent.
5
Real data
We demonstrate Intensity Profile Projection on a dataset containing the face-to-face interactions of the pupils of primary school in Lyon over two days in October 2009 (53). During the study, discreet radio-frequency identification devices were worn by 232 pupils and 10 teachers which recorded their face-to-face interactions. When two participants were in close proximity over interval of 20 seconds, the timestamped interaction event was recorded. The school contains five year groups, each divided into two classes, and each class has assigned room and an assigned teacher. The school day runs from 8:30am to 4:30pm, with lunch break from 12:00pm to 2:00pm, and no data was gathered contacts taking place outside the school during sports activities. For more details about the study and dataset, we refer the reader to [53]
We apply Intensity Profile Projection to the data corresponding to each day of the study using kernel smoother with an Epanechnikov kernel, choosing a bandwidth of 5 minutes and computing 30 dimensional trajectories.
8

---

Day 1
Day 2
&
9:00am
11:00am
1:00pm
3:00pm
5:00pm
9:00am
11:00am
1:00pm
3:00pm
5:00pm
Clases 2
9:00am
11:00am
3:00pm
5:00pm
9:00am
11:00am
1:00pm
3:00pm
Classss &
9:00am
11:00am
1:00pm
3:00pm
5:00pm
9:00am
11:00am
1:00pm
3:00pm
5:00pm
Clasoes & 48
9:00am
11:00am
1:00pm
3:00pm
5:00pm
9:80am
11:00am
1:00pm
3:00pm
5:00pm
classss 51
9:00am
11:00am
1:00pm
3:00pm
5:00pm
9:00am
11:00am
1:00pm
3:00pm
5:00pm
time of day
time of day
class 1A
class 2A
class 3A
class 4A
class 5A
class 1B
class 28
class 3B
class 4B
class 5B
Figure 3: Ome-dimensional PCA visualisation the 30-dimensional node representations for pairs of classes in the same year group. The solid lines show the average trajectory for each class, and the dashed line show one standard deviation above and below.
9:30am, Day
(b) 9:30am, Day
(c) 12:30pm. Day 2
(d) 3:30pm, Day
40
-20
-40
-20
-70
20
-20
-20
20
class
class 2A
class 3A
class 4A
class 1B
class 2B
class 38
class 4B
class 5B
Figure 4: Two-dimensional t-SNE visualisation of the 30-dimensional node representations of all
pupils and teachers evaluated at 9:30am on Day 1, and 9:30am, 12:30pm and 3:30pm on Day 2.

---

To visualise the node trajectories, we first rescale them to have unit norm, which has the effect of removing information about the "activeness" of node from its representation (see, for example, (44). and apply two dimension reduction techniques. The first is principal component analysis (PCA), which we adapt to our dynamic setting by projecting the (centered) representations onto the the direction of maximum average variance over the time domain. This visualisation gives us temporally coherent view on the trajectories (more details are given in Section of the appendix). In Figure 3 we visualise the trajectories of each pair of classes in each year group using PCA, and for clarity, we just plot the average trajectory for each class, along with one standard deviation above and below.
The second t-Distributed Stochastic Neighbor Embeddiig (t-SNE), a popular non-linear dimension- reduction tool which provides enough flexibility to visualise the whole set of representations at each point in time. Figure B throughout the study. In Section B of the appendix, we include analogous figures for aligned spectral embedding and Omnibus embedding for comparison.
shows 1-SNE visualisations of the node representations at collection of times
Figure 3 clearly shows the mixing of classes during the lunch hours, and from Figures 4 we see that the the representations are much more fragmented during the lunch hour (12:30pm, Day 2) than they are during lessons at the other times, where they form tighter clusters corresponding to classes.
While it is reassuring that the geometry of the trajectories reflects the known class and timetable structures of the school, it also allows us to uncover structure in the data that was not known from the report the study. For example, classes 5A and 5B (olive and cyan, respectively) merge into single cluster at approximately 9:30am on Day 1, and classes 3A and 3B (brown and pink, respectively) do the same at approximately 9:30am on Day 2. One might conjecture that this corresponds to a joint lesson, which is taken by the students of both classes in a year group.
6 Discussion
We have presented an algorithmic framework to learn continuous-timg (ow-dimensional trajectories representing the evolving behaviours of nodes in a dynamic network. We view our framework as providing platform on which novel inference procedures can be developed, particularly combining graph and temporal concepts. For example, in dynamic networks with continuously evolving community structure, it might be interesting to develop procedures for detecting branching points (see bifurcating block model example, Section or measures of polarisation and cohesion the network via the velocities of the trajectories. More generally, we believe there is much left to understand and exploit the time-evolving topology and geometry of these representations.
A limitation of our framework is the need for bandwidth and dimension selection. These decisions are difficult because they are trade-offs, bias versus variance in the case of bandwidth selection (as seen here), and statistical versus computational in the case of dimension selection (see e.g. (551). In the presence of a specific supervised downstream task, both decisions could be assisted by cross- validation. In unsupervised settings with ceasonably-sizes networks, our method is very fast, allowing expedient exploration of different choices.
Our theory suggests selecting a dimension which corresponds to an "eigengap" in the spectrum. In practice, it is possible than the spectrum does not decay quickly and no eigengap is present. This likely corresponds to the violation of Assumption Some possible solutions are to take an entrywise transformation of the intensity profiles, such as the square root, to temper heavy tails (56) or to employ robust subspace estimator, such as robust PCA [57].
Our method might be viewed as dynamic analogue of adjacency spectral embedding for static graphs (58) and, as a result, in future research it could be profitable to find dynamic analogues of other variants of spectral embedding, e.g. applying Laplacian normalisation [59.6] or regularisation [54 62].
We believe improved dynamic network analysis can be used for societal good, in applications such as cyber-security, or combating human-trafficking, fraud, and corruption. However, one should also be aware of the risks, particularly to individual privacy and targeted influence.
10

---

References
[1] Mark EJ Newman. Modularity and community strueture in networks. Proceedings of the National Academy of Sciences, 103(23):8577-8888 2006.
[2] Emmanuel Abbe, Community detection and stochastic block models: recent developments. The Journal of Machine Learning Research, 18(1):6446-653, 2017.
[3] Larry Wasserman. Topological data analysis. Annual Review of Statistics and Its Application, 5: 501-532, 2018.
[4] Patrick Rubin-Delanchy. Manifold structure in graph embeddings. Advances in Neural Infor- mation Processing Systems, 33:1 1687-11699, 2020.
[5] Avanti Athreya, Minh Tang, Youngser Park, and Carey E. Priebe. On estimation and inference in latent structure random graphs. Statistical Science, 36(1):68 88, 2021. doi: 10.1214/ 20-STS787. URL BtLL/htps::/Cot :ttttp:://oi..or/11..122222
[6] Patrick O Perry and Patrick Wolfe. Point process modelling for directed interaction networks. Journal of the Royal Statistical Society; Series B: Statistical Methodology, pages 821-849, 2013.
[7] Guotong Xue, Ming Zhong, Jianxin Li, Jia Chen, Chengshuai Zhai, and Ruochen Kong. Dynamic network embedding survey, Neurocomputing,
422.:ppttnnn.4221122222222 472:212-223, 2022.
[8] Paulo Rauber, Alexandre X. Falcão, and Alexandru C. Telea. Visualizing time-dependent data using dynamic t-sne. EuroVis '16 page 73-77, Goslar, DEU, 2016. Eurographics Association.
[9] Ian Gallagher, Andrew Jones, and Patrick Rubin- Delanchy, Spectral embedding for dynamic networks with stability guarantees. Advances in Neural Information Processing Systems, 34: 10158-10170, 2021.
[10] Andrew Jones and Patrick Jubin-Delandhe The multilayer random dot product graph. arXiv preprint arXiv:2007.10455 2020.
[11] Emmanuel Abbe, Jianqing Fan, Kaizheng Wang, and Yiqiao Zhong. Entrywise eigenvector analysis of random matrices with low expected rank. Annals of Statistics, 48(3):1452, 2020.
[12] Lihua Lei. Unified 8 eigenspace perturbation theory for symmetric random matrices. arXiv preprint arXiv: 99990049995 2019.
[13] Jianqing Fan, Weichen Wang, and Yiqiao Zhong. An foo eigenvector perturbation bound and its application to robust covariance estimation. Journal of Machine Learning Research, 18(207):
1 1-42, 2018.
[14] Fangzheng Xie. Entrywise limit theorems of cigenvectors for signal-plus-noise matrix models with weak signals. arXiv preprint arXiv:2106. 099400. 2021.
[15] Jörg Gross and Carsten KW De Dreu. The rise and fall of cooperation through reputation and group polarization. Nature Communications, 10(1):776, 2019.
[16] Joan-Maria Esteban and Debraj Ray. the measurement of polarization. Econometrica: Journal of the Econometric Society, pages 819-851.
[17] Pedro Szekely, Craig A Knoblock, Jason Slepicka, Andrew Philpot, Amandeep Singh, Chengye Yin, Dipsy Kapoor, Prem Natarajan, Daniel Marcu, Kevin Knight, et al. Building and using a knowledge graph combat human trafficking. In The Semantic Web-ISWC 2015: 14th Interna- tional Semantic Web Conference, Bethlehem, PA, USA, October 1 11-15, 2015, Proceedings, Part 14, pages 205-221 Springer, 2015.
[18] Darren Edge. Revealing the hidden structure of corruption. https: / /ww
mi crosoft. cwm/ coom/en-us/reee htap:///////eeerrr////o///aaaaaa hasteaiieeeeeeeeeeeeeeeeeeeee//e htw.lliatenncncees.... articles/
-t -td.ddeesttnnttte-----ooooooppoo 2021. Accessed: 12 May
the-
2023.

---

[19] Carter T Butts. 4. relational event framework for social action. Sociological Methodology, 38 (1):155-200, 2008.
[20] Duy Vu, David Hunter, Padhraic Smyth, and Arthur Asuncion. Continuous-time regression models for longitudinal networks. Advances in Neural Information Processing Systems, 24, 2011.
[21] Francesco Sanna Passino and Nicholas A Heard. Mutually exciting point process graphs for modeling dynamic networks. Journal of Computational and Graphical Statistics, pages 1-15, 2022.
[22] Christopher DuBois, Carter Butts, and Padhraic Smyth. Stochastic blockmodeling of relational event dynamics. In Artificial Intelligence and Statislics, pages 238-246. PMLR, 2013.
[23] Marco Corneli, Pierre Latouche, and Fabrice Rossi. Block modelling in dynamic networks with non-home non-homogeneous poisson processes and exact icl. Social Network Analysis and Mining, 6: 1-14, 2016.
[24] Catherine Matias, Tabea Rebafka, and Fanny Villers. A semiparametric extension of the stochastic block model for longitudinal networks. Biometrika, 10563):665-680, 2018.
[25] Sharmodeep Bhattacharyya and Shirshendu Chatterjee. Spectral clustering for multiple sparse networks: I. arXiv preprint arXiv: 10%554 2018.
[26] Ruthwik Junuthula, Maysam Haghdan, Kevin S Xu, and Vijay Devabhaktuni. The block point process model for continuous-timn event-based dynamic networks. In The World Wide Web Conference, pages 829 839, 2019.
[27] Makan Arastuie, Subhadeep Paul, and Kevin Xu. CHIP: a Hawkes process model for continuous- time networks with scalable and consistent estimation. Advances in Neural Information Pro- cessing Systems, 33:16983-16996, 2020.
[28] Jiasen Yang, Vinayak A Rao, and Jennifer Neville. Decoupling homophily and reciprocity with latent space network models. In UAI, 2017.
[29] Zhipeng Huang, Hadeel Soliman, Subhadeep Paul, and Kevin S Xu. A mutually exciting latent space hawkes process model for continuous-time networks. In Uncertainty Artificial Intelligence, pages 863-873. PMLR, 2022.
[30] Giang Hoang Nguyen, John Boaz Lee, Ryan Rossi, Nesreen K Ahmed, Eunyee Koh, and Sungchul Kim. Continuous-time dynamic network embeddings. In Companion proceedings of the Web Conference 2018, pages 969-976, 2018.
[31] Riccardo Rastelli and Marco Corneli. Continuous latent position models for instantaneous interactions. arXiv preprint 1rXiv:2103.17146, 2021.
[32] Igor Artico and Ernst Wit. Fast inference of latent space dynamics in huge relational event networks. arXiv preprint arXiv:2303. 2023.
[33] Wenjie Fu, Le Song, and Eric Xing. Dynamic mixed membership blockmodel for evolving networks. In Proceedings of the 26th annual International Conference Machine Learning, pages 329-336, 2009.
[34] Kevin S Xu and Alfred O Hero. Dynamic stochastic blockmodels for time-evolving social networks. IEEE Journal of Selected Topics in Signal Processing, 8(4):552-562, 2014.
[35] Catherine Matias and Vincent Miele. Statistical clustering of temporal networks through a dynamic stochastic block model. Journal of the Royal Statistical Society. Series B (Statistical Methodology). 79(4): 1141, 2017.
[36] James D Wilson, Nathaniel T Stevens, and William H Woodall. Modeling and detecting change temporal networks via the degree corrected stochastic block model. Quality and Reliability Engineering International, 35t 1363, 2019.
12

---

[37] Purnamrita Sarkar and Andrew W Moore. Dynamic social network analysis using latent space models, Advances in Neural Information Processing Systems, 18:1145, 2006.
[38]
Daniel K Sewell and Yuguo Chen. Latent space models for dynamic networks. Journal of the American Statistical Association, 1 (166512)1644-111 2015.
[39] Nial Friel, Riccardo Rastelli, Jason Wyse, and Adrian E Raftery. Interlocking directorates in Irish companies using a latent space model for bipartite networks. Proceedings of the National Academy of Sciences, 1 113(24):6629-666) 2016.
[40] Keith Levin, Avanti Athreya, Minh Tang, Vince Lyzinski, Youngser Park, and Carey E Priebe. central limit theorem for an omnibus embedding of multiple random graphs and implications for multiscale network inference. arXiv preprint arXiv: / 705505 2017.
[41] Fuchen Liu, David Choi, Lu Xie, and Kathryn Roeder. Global spectral clustering in dynamic networks. Proceedings of the National Academy of Sciences, (155)9927-932, 2018.
[42] Joshua Cape. Spectral analysis of networks with latent space dynamics and signs. Stat, 10(1): e381, 2021.
[43] Lun Du, Yun Wang, Guojie Song, Zhicong Lu, and Junshan Wang. Dynamic network embedding: An extended approach for skip-gram based network embedding. In IJCAI, volume 2018, pages 2096-2092, 2018.
[44] Sedigheh Mahdavi, Shima Khoshraftar, and Aijun An. dynnode2vec: Scalable dynamic network embedding. In 2018 IEEE International Conference on Big Data (Big Data), pages 3762-3765 IEEE, 2018
[45] Carl Eckart and Gale Young. The approximation of one matrix by another of lower rank. Psychometrika, 133:2111-118, 1936.
[46] Richard Bruno Lehoucq. Analysis and implementation of an implicitly restarted Arnoldi iteration. Rice University, 1995.
[47] James Baglama and Lothar Reichel. Augmented implicitly restarted Lanczos bidiagonalization methods. SIAM Journal on Scientific Computing, 27(1):19-42, 2005.
[48] Keith Levin, Fred Roosta, Michael Mahoney, and Carey Priebe. Out-of-sample extension of graph adjacency spectral embedding. In International Conference Machine Learning, pages 2975-2984 PMLR, 2018
[49] Emmanuel Candes and Benjamin Recht. Exact matrix completion via convex optimization. Communications of the ACM, 55(6):111-11,, 2012.
[50] Avanti Athreya, Donniell E Fishkind, Minh Tang, Carey E Priebe, Youngser Park, Joshua T Vogelstein, Keith Levin, Vince Lyzinski, and Yichen Qin. Statistical inference on random dot product graphs: a survey. The Journal of Machine Learning Research, 18(1):8333-848) 2017.
[51] Emmanuel Abbe and Colin Sandon. Community detection in general stochastic block models: Fundamental limits and efficient algorithms for recovery. In 2015 IEEE 56th Annual Symposium on Foundations of Computer Science, pages 670-688. IEEE, 2015.
[52] Peter W MacDonald, Elizaveta Levina, and Ji Zhu. Latent space models for multiplex networks with shared structure. Biometrika, 109(3):683-766, 2022.
[53] Juliette Stehlé, Nicolas Voirin, Alain Barrat, Ciro Cattuto, Lorenzo Isella, Jran-François Pinton, Marco Quaggiotto, Wouter Van den Broeck, Corinne Régis, Bruno Lina, et al. High-resolution measurements of face-to-face contact patterns in a primary school. PloS one, 6(8):e2317,
2011.
[54] Tai Qin and Karl Rohe. Regularized spectral clustering under the degree-corrected stochastic blockmodel. Advances in Neural Information Processing Systems, 26, 2013.
[55] Matthias I Löffler, Anderson Zhang, and Harrison H Zhou. Optimality of spectral clustering in the gaussian mixture model. The Annals of Statistics, 49(5):2506-2534) 2021.
13

---

[56] Ian Gallagher, Andrew Jones. Anna Bertiger, Carey Priebe, and Patrick Seeinn-ellnncy.... Spectral embedding of weighted graphs. Journal of the American Statistical Association, (justacceepedd): 1-19, 2023.
[57] Emmanuel J Candès, Xiaodong Li, Yi Ma, and John Wright. Robust principal component analysis? Journa of the ACM (JACM), 38,
[58] Daniel Sussman, Minh Tang, Donniell Fishkind, and Carey Priebe. A consistent adjacency spectral embedding for stochastic blockmodel graphs. Journal of the American Statistical Association, 107(499):1119-11 2012.
[59] Fan RK Chung. Spectral graph theory, volume 92, American Mathematical Soc., 1997.
[60] Ulrike Von Luxburg. A tutorial on spectral clustering. Statistics and Computing, 2007.
[61] Minh Tang and Carey Priebe. Limit theorems for eigenvectors of the normalized laplacian for random graphs. The Annals of Statistics,, (6(5):2360-24111 2018.
[62] Yilin Zhang and Karl Rohe. Understanding regularized spectral clustering via graph conductance. Advances in Neural Information Processing Systems, 31, 2018.
[63] Francesco Sanna Passino, Nicholas A Heard, and Patrick Rubin- Jubin-Delanchy. Spectral clustering on spherical coordinates under the degree-corrected stochastic blockmodel. Technometrics, 64 (3):346-357, 2022.
[64] Laurens Van der Maaten and Geoffrey Hinton. Visualizing data using t-sne. Journal of Machine Learning Research, 9(11 1), 2008
[65] George Linderman, Manas Rachh, Jeremy Hoskins, Stefan Steinerberger, and Yuval Kluger. Fast interpolattonn based t-sne for improved visualization of single-cell ma-seq data. Nature Methods, 16(3):243-245. 2019.
[66] Fangzheng Xie. Entrywise limit theorems of eigenvectors and their one-step refinement for sparse random graphs. arXiv preprint arXiv:2 200. 09840, 2021.
[67] Stanislav Minsker. On some extensions of bernstein's 's inequality for self-adjoint operators. Statistics & Probability Letters, 127:111-119 2017.
[68] Afonso S Bandeira and Ramon Van Handel. Sharp nonasymptotic bounds on the norm of random matrices with independent entries. Annals of Prohability, 44(4):2479-2500 2016.
[69] Yuxin Chen, Yuejic Chi, Jianqing Fan, Cong Ma, et al. Spectral methods for data science: A statistical perspective. Foundations and Trends@ in Machine Learning, 14(5):566 806, 2021.
[70] David Pollard. Empirical processes: Theory and applications. In NSF-CBMS Regional Confer- ence Series in Probability and Statistics, pages i-86. JSTOR, 1990.
14

---

Appendix to "Intensity Profile Projection: A Framework for Continuous-Time Representation Learning for Dynamic Networks"
A Details of the simulated example and method comparison of Sections 4.1 and 4.2
We simulate data according to the following generative process, which might be viewed as describing as a dynamic, two-community stochastic block model. Assign to each node a variable 2; {1,2] denoting its community (which does not change with time). If nodes i and are in the same community, i.e., 23 = 2j. the point process Eig follows a homogeneous Poisson process with (fixed) intensity Nj-- Otherwise, Eij follows inhomogeneous Poisson process with intensity
no exp(m - (1))
< 811
X1y(1) =
81 S t $2
no exoepp-mic -t)
$2)}
82,
where 0 < $1 $2 < T. This model describes two communities gradually coming together until fully merging by time $1, before splitting at time $2 and then gradually drifting apart. We simulate from this model using the parameters T 1, n 100, 21,..., 250 1, 251 2100 2, 170 = 100,
?/ = 10, $1 = 0.3 and 82 0.7.
In our method comparison, we used an embedding dimension of d 2 for all methods, unless otherwise stated. For the discrete-time methods, we construct a series of snapshots of the continuous-time network, each a weighted static network whose edge weights are the number of events which occur on edge in the corresponding time window. The selection of hyperparameters for each method is outlined below:
Intensity Profile Projection (histogram): We used bin size of
1 M = 20
Intensity Profile Projection (kernel smoothing): We used a Epanechnikov kernel with bandwidth 0.1 and applied the approximate Intensity Profile Projection algorithm with B = 20. Different values of bandwidth gave similar results in terms of embedding structure; we chose this bandwidth achieve the desired smoothness.
CLPM 31 The dimension is automatically computed by the algorithm d = 2. The hyperparameters are chosen equal to the ones used in "Simulation C" in [31 which is similar simulated example with two communities. We used 19 changes point which correspond to 20 windows. The implementation was obtained from the Github repository http:////ggtuuub..ooo/arrrrgggge cwt.///aaooogne///.....
The dimension is automatically selected by the algorithm as d 2. The smoothing parameter is chosen with cross-validation which results in equivalent log- likelihood values for from 0.00001 to 0.001. We choose 0.001 which is the larger value for which the algorithm converges. The implementation was obtained from the Github repository https //ggtttttuuug com/ www.nn/Pisscaa.
PisCES 41):
The dimension is automatically selected by the algorithm as d = 2 for all windows except windows to 16 for which we set missing the second dimension to zeros. The implementation is obtained from the multinest R package (available on CRAN) and hyperparameters are set to their default values.
MultiNeSS 52
selected. For these windows,
15

---

9:30am, Day
(b) 12:30pm, Day
3:30pm, Day
2.0
1.5
(d) 9:30am, Day
(e) 12:30pm, Day
(0 3:30pm, Day 2
2.5
2.0
1.0
0.5
1.0
1.5
0.5
1.0
1.5
class
class
class
class
Figure 5: The first two dimensions of the spherical coordinates of the coordinates X,(() using the histogram intensity estimator for times corresponding to the morning, lunchtime and afternoon across both days. The colours indicate classes with black points representing teachers.
(a) 9:30am, Day
(b) 12:30pm. Day
(c) 3:30pm, Day
2.0
0.5
9:30am, Day
(e) 12:30pm, Day
(n) 3:30pm, Day
2.0
1.5
0.5
1.5
2.0
2.5
3.0
1.5
2.0
2.5
3.0
1.5
2.0
2.5
class 1A
class 2A
class 3A
class 4A
class 5A
class 1B
class 2B
class 3B
class 4B
class 5B
Figure 6: The first two dimensions the spherical coordinates of the trajectories (t) using the Epanechnikov kernel smoother for times corresponding to the morning, lunchtime and afternoon across both days. The colours indicate classes with black points representing teachers.
16

---

Day
Day
Cluoses
14 6
3:30am
11 11:0
1:00pm
3:00pm
5:00pm
9:00am
11:00am
1:00pm
3:00pm
5:00pm
Classes
24 2..28
9:00am
11:00am
1:00pm
3:00pm
5:00pmm
9:00aam
11:00am
1:00pm
3:00pm
5:00pm
Classes
30 9:38
9:00am
11:0
1:00pm
3:00pm
5:00pm
5:00am
11:00am
1:00pm
3:00pm
46 &r
8:00am
11 10:00ammm
1:00pm
2:00pm
5:00pm
9:00am
11:00am
1:00pm
2:00pm
5:00pm
Classun 54 $6
9:00am
11 11:00am
1:00pm
3:00pm
5:00pm
9:00am
11:00am
1:00pm time of day
3:00pm
5:00pm
time of day
class
class
class
dlass 4A
class SA
class
class
class
dlass
class
Figure 7: One-dimensional PCA visualisation the 30-dimensional node representations for pairs of classes in the same year group. The solid lines show the average trajectory for each class, and the dashed line shows one standard deviation above and below.
(a) 9:30am, Day
(b) 9:30am, Day 2
(c) 12:30pm, Day 2
(d) 3:30pm, Day 2
-20
-40
-20
-20
-20
class 1A
class 2A
class 3A
class 4A
class 5A
class 1B
class 2B
class 3B
class 4B
class 5B
Figure 8: Two-dimensional t-SNE visualisation of the 30-dimensional node representations of all
pupils and teachers evaluated at 9:30am on Day 1, and 9:30am, 12:30pm and 3:30pm Day 2.
17

---

Day
Day 2
Clocoess
9:00am
11:00am
1:00pm
3:00pm
5:00pm
9:00am
11:00am
1:00pm
3:00pm
5:00pm
Clocses
2A 28
9:00am
11:00am
1:00pm
3:00pm
5:00pm
9:80am
11:00am
1:00pm
3:00pm
5:00pm
Classes 3A 30
9:00am
11:00am
1:00pm
3:00pm
5:00pm
9:00am
11:00am
1:00pm
3:00pm
pm00
Clacsrss
6 46
9:00am
11:00am
1:00pm
3:00pm
5:00pm
9:00am
11:00am
1:00pm
3: 3:00pm
5 5:00pm
Classes 5B
9:00am
11:00am
1:00pm
3:00pm
5:00pm
9:00am
11:00am
1:00pm
3:00pm
5:00pm
day
day
3A
class
class
Figure One-dimensional PCA visualisation the 30-dimensional node representations for pairs of classes in the same year group, obtained using aligned spectral embedding. The solid lines show the average trajectory for each class, and the dashed line shows one standard deviation above and below.
18

---

Day 1
Day 2
1A& 18
9:00am
11:00am
1:00pm
3:00pm
5:00pm
9:00am
11:00am
1:00pm
3:00pm
5:00pm
Closes
&
9:00am
11:00am
1:00pm
3:00pm
5:00pm
9:00am
11:00am
1:00pm
3:00pm
5:00pm
Closses
3
9:00am
11:00am
1:00pm
3.00pm
5:00pm
9:00am
11:00am
1:00pm
3.00pm
Clacees 4A & 48
1:00pm
3:00pm
5:00pm
9:00am
11:00am
1:00pm
3:00pm
Classes & 58
9:00am
11:00am
1:00pm time of day
3:00pm
5:00pm
9:00am
11:00am
1:00pm time of day
3:00pm
5:00pm
class 1A
class 2A
class 3A
class 4A
class 5A
class 1B
class 28
class 3B
class 4B
class 5B
Figure 10: One-dimensional PCA visualisation of the 30-dimensional node representations for pairs of classes in the same year group, obtained using the Omnibus spectral embedding. The solid lines show the average trajectory for each class, and the dashed line shows one standard deviation above and below.
19

---

B Additional real data analysis
We provide further experiments and details of the Intensity Profile Projection analysis of the Lyon primary school dataset described in Section [5] As a comparison to the analysis in the paper, we apply Intensity Profile Projection to the data corresponding to each day of the study with a histogram intensity estimator, choosing a bin size of minutes and computing 30-dimensional trajectories.
Figures 5 6 the first two spherical coordinates (33 of the trajectories obtained using the histogram intensity estimator and the Epanechnikov kemel smoother, respectively. The six plots correspond to the morning, lunchtime and afternoon across both days.
Figure 7 (equivalent to Figure visualises the trajectories of each pair of classes in each year group using PCA where we plot the average trajectory for each class, along with one standard deviation above and below. Since every trajectory using the histogram intensity estimator is piece- wise constant, so are the resulting PCA averages. The pairs of trajectories merge and split in a similar way to those obtained using the kernel smoother.
Figure shows t-SNE visualisations of the node representations at a collection of times throughout the study. The plots are very similar to the equivalent Figure for the kernel smoother with almost identical clusters of students before and after lunch, albeit placed differently by the t-SNE algorithm.
We apply the aligned spectral embedding method and the Omnibus spectral embedding method, described in Section 4.2 on the main text, to snapshots of the data. Figures and TO visualises the respective trajectories of each pair of classes in each year group using PCA where we plot the average trajectory for each class, along with one standard deviation above and below. These representations fail to capture the behaviour of the network and fail in the same ways as they did for the simulated data in Section 4
C Visualisation
In this section, we give short overview of the two dimension reduction techniques employed for visualisation in this paper.
For the trajectory visualisation in Figures 2 and we use a principal component analysis which we extend to the dynamic setting by computing projection using the leading eigenvectors of the average covariance matrix, which we apply to the (globally centered) trajectories. This has a similar flavour to our Intensity Profile Projection algorithm, and since we reduce dimension using a common projection, it gives a temporally coherent view of the trajectories.
The second visualisation technique we apply is t-SNE (64), using the FIt-SNE implementation [65], which we used to obtain Figure 4 This visualisation method is not naturally extended to dynamic data, so we initialise the algorithm using the aforementioned dynamic extension PCA, which results in the visualisations at different times being approximately aligned.
D Computational complexity
In this section, we given brief discussion of the computational complexity of Algorithm [ Suppose we use a kernel with finite support of width 2/, then the expected time complexity of the kernel intensity estimate for single edge at a single point in time is O(hAmax). Evaluation of A(t) is hhAha and of the matrix (3) is Q(Bn2h. h-maxx) Consider the setting of Corollary | where (t)(t) = p for all i, and L x PLa where μ4, d, Lo x 1 are fixed, our theory suggests choosing h (np) -1/3, and if we suppose np log" n, the sparsest regime our theory allows, then evaluation of 3 is O(Bn log" m)), i.c. lin- in n.
In practice the top singular vectors of (3) can be computed using the Augmented Implicitly Restarted Lanczos Bidiagonalization algorithm [47] implemented in the irlba package in R, or the irlbpy in Python. The time complexity of this algorithm has not been studied theoretically although in practice it can be incredibly fast. For example, the package author performs a simulated experiment in which they compute the first 2 singular vectors of sparse 10M x 10M matrix with IM non-zero entries which takes approximately 6 seconds on computer with two Intel Xeon CPU E5-2650
20

---

processors (16 physical CPU cores) running at 2 GHz equipped with 128 GB of ECC DDR3 RAM (see htwptt::/bbwwwwss.ggttettee githuu...gtt..........g. hot.://aaooooooooaaaaaaaaaaaaaaa html).
E
Proof of Lemma
We begin by writing
[ (1; V)d.
arg min VEO(n.d)
(V) arg min VEO(n.d)
VV ( - -A(t dt
arg min VEO(n.d) 0
1=1
(I- VV (V-VV)) dt,
arg min VEO(n.d)
and since (I - VVT) is the projection onto the orthogonal complement of the columns space of V, we have that
| = / (aa
Therefore, minimising R2 (V) is equivalent maximising
[ Σ
Σ V u
| V
where the equality holds due to the invariance of the Euclidean norm under orthogonal transformations. As result, we have
arg min VEO(n,d)
(V( arg max VeO(n,d) 0 1=1
Σ
v A
|
arg max VEO(n,d)
tr (VTVV dt
arg max VEO(n.d)
VT ( (t(() VIEV
arg tr VEO(n.d)
arg tr VEO(n.d)
where the final equality follows from the Courant- Fisher min-max theorem. This concludes the proof.
F Proof of Theorem 1
F1 Prerequisites
F.1.1 Additional notation
In this proof, we use the notation S bn to mean an, S by with overwhelming probability.
F.1.2 Symmetric dilation with change of basis "trick"
Symmetric dilation is proof technique which allows statements about the eigenvectors of symmetric matrix to be easily extended to hold for the singular vectors of (potentially rectangular) asymmetric
21

---

matrix. Let M be an x matrix with non-zerr singular values {or and corresponding orthonormal left singular vectors {urli=1 and right singular vectors f0r Its symmetric dilation is the n X n matrix (with + 122) constructed as
D(M)
M
0
One easily verify that (M has cigenvalues {= and eigenvectors {(u, ) Fi=1 We stack the first d left and right singular vectors into matrices U E I [Raxxd xd and V E RRxxd and stack the first 2d eigenvectors of D(M) into matrix
1
U V
U -V
U
V2
We then have
2 V V
1,,x 2,00 =
| 2,00
and
| = ||
While this standard construction very useful when n1 = f2, it can lead to suboptimal bounds when n2 > n1, or >> n2, due to issue about incoherence, which was first raised in [13. The incoherence of a subspace Uo spanned by the orthonormal columns of a matrix Up E R
1 (U))
|oll2.00
To obtain a good entrywise eigenvector bound under sgnaal-puu-noose matrix model is typically necessary that (U) = 1. Observe that
(U)
11 +
T 11,,,00
111 122
(1U112,0
+ 112
(U)++
11
112
1 (V).
2d
2d
2n1
2n
If (U), (V) x and n1 x n2, then μ (UU) x 1 and is typically possible to obtain good bounds. However, when n2 > n1, we have (() 1, and good bound can typically not be obtained. The imbalance of and ny can cause similar issues when obtaining spectral norm bounds.
This issue be overcome by changing
basis which balances the contribution from its first nı and second n2 elements of each column. Specifically, let лh /2m1/ +22) and T2
/ /2n2/ + n2), and consider the basis E....., En1 +na+ such that
T if E (l....,.1) m13
=
T2ē;
if E {n1 + I,..., n1 + 122,
where {e:)_=t are the standard basis vectors in R. Let III-I1,, denote norm with respect to the column basis [2 121422 then one can verify that
[[ [M]], and
1111, = |2,00 2,00 T V2lx2
As result, if (UU)) / na/d | Un 2.00 then
A(Ú)
V(().
regardless of the relative sizes of n, and 122. We use this symmetric dilation with change-of- "trick" to apply some existing theorems for symmetric matrices to our setting.
F.1.3 Concentration inequalities
In this section, we state a collection of lemmas which will make use of throughout the proof. We begin with tail bound for a Poisson random variable.
Lemma 3. Let X
Poisson(A) A). Then
2 2(A + t/3)
P P(X-- A 2 2 exp
For 2 X,
P[[X 1/
2e
31/8
22

---

The bound can be established by approximating the Poisson distribution with mean À as the sum of l: Bernoulli random variables with mean A/k. applying Bernstein' 's inequality, and taking
Our next result is a concentration bound which adapts Lemma A. of (66] and can be proved using a vector version of the Bernstein inequality (Corollary 4.1 in (67)
Lemma 4. Let X; ~ Poisson( Ac) independently for all i = 1..... n, and suppose Q E Rexd is deterministic matrix whose rows denote Qi- Let Amax == maxie[n] di then with probability
1 28n
Ex, -x,
S 3log* n 11QQ12,,00 + 6 log [Q]|
2
Next, we state a concentration bound for the spectral norm of random matrices with independent entries which appears as Corollary 3.12 in [68] The original statement of this lemma is symmetric random matrices, although we general it to arbitrary random matrices using the symmetric dilation with change-of- of-age----bbasss trick described in Section F.1.2
Lemma 5 (Corollary 3.12 of [68]). Let X be an ny x matrix whose entries Tij are independent random variables which obey
E (xay) = O, and Digl B,
i [m1],j E [m22].
Then there exists universal constant > 0 such that for any 2
t²
P X| 2 4/D S n nxxp
cB2
where
= max
max EE
(
max
Σ
E (zft)
[ma]
==1
j=1
and Fy. 221x/(n1 + n2)-
F.1.4 Weyl's inequality and Wedin's 's sinO theorem
The next two lemmas are classical results matrix perturbation theory. Weyl's inequality shows that the singular values of matrix are stable with respect to small perturbations.
Lemma 6 (Weyl's inequality). Let M, E be n, x f2 real-valued matrices. Then for every 1 sis (nu n2), the ith largest singular value of M and M obey
l (M E) σg (M)l
[
One way to measure the distance between two subspaces 24 and E is via principal angles. Let U, Ü be matrices whose orthonormal columns span 2I and úi respectively, and let (5)) denote the singular values of U U. Then the principal angles {0. [ between IA and i are defined by & cos(O.). Let sin (U, Ū) == diag(sin o,,..., 0a). Another way to measure the distance between and ú is via the difference between the projection operators UUT and ÜÚT, and in fact, these two characterisations are equivalent. Specifically,
sin (U, [ = -
We will use this equivalence without mention throughout the proof. Wedin's 's sinO theorem shows that the singular vectors matrix are stable with respect to small perturbations.
Lemma 7. Let M and M = M + E be two n. x n2 real-valued matrices, and denote by U, Ô (respectively V. V) the matrices whose columns contain d orthonormal left (respectively, right) singular vectors, corresponding the d largest singular values of M and M. Let d aM ddt1 (M) and suppose that |E]] (1 1/V2)6, then
sin + (u. sin e (v, v| <
|
EV
2 E
8
See 669 for proof.
23

---

F.2 Implications of Assumptions 1
We state here some inequalities involving the parameters of our problem which follow from As- sumptions 14 and elementary linear algebra. We will use these facts throughout the proof without mention.
V Vndmax S S nAmaxi
(6)
log 2 KnAmaxi
(7)
logn.
(8)
The inequality (6) holds since S o 1/2 0//2(2)) 2) S vn
Σ II S ndmax and n2mmx log 2 Vndmax
kn.mmax log(ō/ VTXmaz)
ndmax log n
where we invoked Assumption
holds by noting that the previous bound implies
logn and invoking Assumption 5 log n nAmax S log n.
log(6/ (o(((//nnnun VnAmax)
8 follows from 66) since X
F.3 Setup
We begin by defining M equally spaced bins in (0, 1}.
2
M -
Bi 0, M B2 M M
.B
and define the piecewise approximation of Az(1),
(,(1) = M J A,(()d,
LE Bm, m E [M].
Bm
We then define f1. .......... tm (0, 1] such that (t)(t) Xt for all Bm, which exist by the continuity of (t), and define the piecewise constant approximation of Yt(t) as Y,(t) UT (t).
Our strategy to obtain the bound in Theorem is to decompose into bias and variance terms:
W P.(1) - x:00112 max sup W w. 81(1 - 7:(0 + max sup W2 Y2(1) Y:01l2 (JT
max sup (JE(M) IET
(ET
variance
bias
Section F66 is dedicated to bounding the bias term, and the rest of this section is dedicated bounding the variance term. Define the unfolding matrices and A (without arguments) and their (thin) singular value decompositions
Ä (A(t)) A(tas) USVT + +0.SVV.
Á = (A(t1) A(tM)) = USV + +Ú_S_VI.
Then one has that for / E Bm, m [M],
(t) A(m Ū Sn Y()) = A(t Ü= V..S
where V m. V.. m denote the with blocks of V and V respectively. Therefore it follows that,
W2 T.00 VSW. - VS -S
max sup CJEIMi (ET
2.00
For ease of exposition, drop the subscript W, in this section. Our bound is based on the following decomposition of VS VSW.
Proposition 1. We have the decomposition
VS VSW V(V VS-SW) SW)
(9)
+(I VI-VVJAT ÚW)
(10)
+ +[- (Ā U)
(11)
+ + ( (
(12)
24

---

Proof of Proposition 7 We begin by adding and subtracting terms to obtain
VS VSW = VS - VV V(V VS SW).
Then, noting that VS Ā and (I VV 0, we have
VS - VV
V$ = - VV A
= -1- V
=(1 (I- VV VV-VV)(A -A-AA - (I-VV'A'O
= (1 VV-VV - ( UW)
11-0-(1-V''A'
Next, we decompose (I- VVT) T by adding and subtracting terms to obtain
VVVVVV - -
(I
UW
VV (I- A) (Ū UW)
12
E4 Technical propositions
We now outline series of technical propositions which require bound terms - which prove in Section G
Our first proposition a 1-norm and spectral norm bound for A.
Proposition The bounds
A A- S M log n.
A A S V Mnnmmxx
hold with overwhelming probability.
The spectral norm bound is obtained using Lemma 5 and the 1-norm bound obtained via an application of the classical Bernstein inequality. The next proposition provides control the singular values of A.
Proposition 3. Let oz(.) denote the ith ordered singular value of a matrix. The singular values of Å satisfy
Moa(EE) S Sa((A) S ai (A) S (M)
The result is obtained using Weyl's inequality. The next proposition provides control of the spectral norm of Q QTA A)R, where Q R are conformable, deterministic unit-norm matrices. Proposition 4. For conformable, deterministic unit-norm matrices Q, R. the bound
(AT(A-A)R]], QQ (Ā (a x)R||, S M log
(13)
holds with overwhelming probability.
The proof of Proposition 4 emmloys classical e-net argument to the spectral norm of appropriately constructed symmetric dilation matrix.
The next proposition states that both the matrices Ü T Ü and V V are well approximated by common orthogonal matrix.
Proposition 5. There exists orthogonal matrix W such that
|
VTdmmax
V
W <
Vndmax
hold with overwhelming probability.
25

---

To prove Proposition we empoy the Wedin sine theorem to obtain bound on | W We then obtain bound | U
V V 22, and combine these bounds to establish the proposition.
The next technical tool require is the ability to "swap" W, W.S and S. Proposition 6. The bound
WS sw/, s M log
holds with overwhelming probability.
This result follows by applying the previous propositions to an appropriately constructed decomposi- tion.
Part of the challenge of obtaining good bound the term 12) is that (A Á) and (Ü UW) are dependent, and this dependence must be decoupled order to apply the standard suite of matrix perturbation tools. For 1, ...., let
N (i.j) : i m or E {m + (e 1)n, E [M]}
and construct the auxiliary matrices A (2),
A (m) defined by
(m)
Ai
if (i,j) & N
(14)
Aij
if (i,j) E Nm.
In words. ((r)) the matrix obtained by replacing the mth row and columns each of its blocks with its expectation. In this way, the mth row of (Å - A) and (m)) independent. Let (mm)
denote the matrix of leading left singular values of A (m)
We apply result due (mm) ((m) - U112,00-
which provides &2,00 control of | ( 12,,00 and
Proposition 7. The bounds
| 2,00
|
(i(m) - -
Admax v dog
12,00 2,00
2,00
hold with overwhelming probabilit.
In addition, require control the spectral norm difference between the projection matrices
((m) (Ū(m))T and the projection matrices ÜÜT and ÜÜT, which is provided in the following proposition.
Proposition 8. The bounds
(mm) ((())T ÜÜT
A1Amax
(15)
2
[Amax d 3/2
3/2
(0m)) (m(m)) T
(16)
82
hold with overwhelming probability.
The proof of Proposition g requires delicate "leave-one-out" argument.
F.5
Bounding terms 9 12
Firstly observe that
VndAmaxd
|| = AATSS-I '12,00 S /
| 2,00
| |- <
(M)
26

---

and therefore term can be bounded as
V(V VS SW)
V12,,00 ( - w S + /w sw(,)
-W
2,00
V2 AV
Vnddmaocf Vndmax (Mo + M 1ogg3/2
V /Moa(2)
nV Mddmmaxf UI
X HVV M Mamax
where the third inequality follows from Assumption that VM log n S nAmax. and the /o)(E)) and the fourth inequality follows from Assumption 3 that
definition (0(l///v/nna vndmaxx Z knAman.
ö log n 2 ö log(6/
To bound 10). we first apply Wedin's sin theorem to obtain
A A
- uw
|| sin @ 00|l2 S
@
V2
Mndmax
Vndmax
ad(A) (A)
VM5
Then, we use Assumption 2 to obtain the bound
|| 2,00 A (I (1) S max suprspp S selnl LET
pd log 5/2 n.
log
Putting these two bounds together, we bound (10) as
O ÜW)
s |a VV-VVTJM J
1
- - UW
Ú-ÜW
PV ddmax 302
2,00
2,00
To bound term (I),, we set Ā
and note that each column of M --EE contains independent Poisson random variables with means no greater that M M- We will use Lemma 4 to bound the rows EÜ as
E U = Σ en S M log"
2.00 + MM Amax logn Ü n 0
12.00
j=1
M log² MMndmax log
| 2,00
S v M Amaxii log n
Ü
2,00
MV Momax
where the third inequality uses Assumption 4 and union bound over i = I,...,n, Therefore, we have
(Ā - A)
2,00 S M do n.
(17)
Noting that II- VV
I VVT
112 S bound 11
(I-
I
v
(A-A)
(Λ-A)
U
MV MAnaed log
2,00
2,00
Finally, we bound term 12 Let A A (m) denote the auxiliary matrices described in (4).
and (mm) denote the matrix of leading left singular values of Am)) We can then decompose the
Euclidean norm of (A A) (Ū UW) as
A A) (Ō- ÜW) /| s
A A) Ü(W i OTUilL
(18)
+ - A) U -
(19)
(~ U)
+
(20)
27

---

The first term
18 is bounded as
A A) Ú( - U ÚTÚ)
S A ~
w ū
2,00
udmax v dn n
nnmmx
M MnAmar log n
/M max UV dlog
X
3/2 n
82
2 /vV M Amaxd log 3/2 n.
To bound the second term 19). we employ Proposition g to obtain
a - A) (0Ú U- ( (m)
< A A
00 - -O(m) (mm))
M.maxx v dn log 3/2
MTidmxx
82
V Mu
vdn 3/2 log" 3/2
n
52
S UV M Amaxd 1/2 n
3/2
We now set about bounding the third term 20 Let 22 E.n. denote singular value decomposi- tion of (m))) T Ü, and set ((m) == S2 Let mlm) denote the principal angles between the
column spaces of (mm) and U defined by <!m)) (Üm)) Ú. We invoke Wedin's 's theorem to show that
cos()!")) where {. are the singular values of
(()) - (ū(m) U
| B11 = max - mm) = max (1 cos m
=
2
EE|d]
sE|d]
^ (m) A
S max(1 - cos? p(m))
max sin 2 mfm) xE[d]
Mndmax
ndmax
iE(d]
(a. J2+1 (A))
M8 8²
We define (m)) = (m)) (Ū (m)) TT UUT and note that H (m)) is independent of (A A mm. and that
H m)
[T(m) W((m) - T// 2,00
+
(mm)
(Ū(m))TT W(m))
2,00
2,00
2 24
( w(()) 0
+
[m))
2,00
2,00
Wmmax vdn logn
Then, using Lemma 4 we have that
|(A A) H m))
M log² n H())
+ Amaxx log n
H (m)
2,00
M log n1maax
H (m)
+
Amaxd logn
H
2,00
udmax v dn log n
Amaxel log ndmax
M log n.
8
8
v M Munii
3/2
3/2og3/2
3/2
8
Md
log
n.
28

---

Combining these bounds and taking union bound over m E [nx], we have
A ()-
S PV MdAmax log 5/2 n2,
2,50
and the term
12) is bounded as
[I- O -
I VV (A- A)
-
AA
2,00
2,00
2/
MV M Ammaxe log
Combining the bounds on 9 12
we have
W W2 7100
VSWT - VS
PV Mdd log n,
max sup
(JE(m] SE
2,,00
which completes the proof.
F.6
Controlling the bias term
F.6.1 Edge level bias
We begin by studying the edge-level bias of the histogram intensity estimator. Let Piy(t)
So )ty(8 (8) ds denote the cumulative intensity of edge j. Now we have, for t E Br,
(t) M / (s)(s) dt
( ((-1)}
M
Pij M - Pig
M
((x))
Pij
M
Pij
M
M
(!))
for some t" E Bf, which follows by an application of the mean value theorem, where d'(t) = (t)(t).
We then apply the L-Lipschitz continuity of X13 (t) obtain
Aty(t) (t)( = [ig(t*) Aeg (t)
$ L lt" t]
M
F.6.2 A subspace perturbation bound
Define the operator A: (T + 12a) 1 by
Ac(-) A(1)(1)) de
and define the operator :R" + (T + RR)) by
Aru A(-)u.
Then Σ AA" since
A AA'u=A(A()u) / A? dt Eu.
Denote its eigenvalues of.... 02 and its corresponding orthonormal eigenvectors ui,..,, Uyy and define v²(-) = A(-)u,/E, for all 1,...,n. n. Then, A-)) admits the (functional) singular value decomposition
m
A(- Σ o,
29

---

Define A and its corresponding parameters analogously. By definition,
nL
A - sup SET
(t) (t) S sup max Aty(t) - Xcy(t) - LE7 JJEjn]
M
Therefore, by (a functional version of) Wedin's 's sin θ theorem
A All nL
|UW U11 S
0d Ja+1
M6
F.6.3 Controlling the bias term
Combining the above bounds, we have that uniformly for all i, t.
[(()), 2221)22 = ||()) A())
|(0))2, | -Ul, ||A(t) A ||1|
n Amax
vnL
+
M6
M
n/ X
M6
where the final inequality follows from the fact that 8 S n2mmax-
G Proofs of the technical propositions
G.1 Proof of Proposition 2
We have that -1 times the lower-triangular elements of each block of Á are independent Poisson random variables with mean given by M- times the lower-triangular elements of each block of A. Define the matrices L and A with the upper and lower triangles, respectively, of each block set zero, and the diagonals of each block halved, and define AL and A similarly, that ÀL (respectively - U has independent Poisson entries with means M-1 A (respectively M- and A A = A - - + (A" ).
We condition on the event that (Å A)) S log for all i,j. which occurs with overwhelming probability by Lemma and union bound. Now, we employ Lemma 5 with B = M log and = Mndmax to obtain
P (( AL 2
+t) s nexp (-
t²
c(M log n)
Setting t = M log 3/2 1/2 n. we have that
AL A S Mndmaax M log S / /Mndmax
where the final inequality follows from Assumption A A and combine the with the triangle inequality:
We obtain an analogous bound for
^ A s A
1411
+
AU
A
Mn
2
We now establish bound on A
We condition on the event A A S M log n. for all i.j. which occurs with overwhelming probability due to Lemma and union bound, and note that
we have Σn= B(A )2 < Mn.mux. Then, by the classical Bernstein inequality, we have for
any >
1
P
A
Ayo
2t < 2exp
2 (M + LM log n/3)
30

---

and setting t = vnM Amax log n. we obtain
Σ
A Azs/ S nMM Xmax log n.
==1
union bound establishes that
P 11, 2 V Amax log nt
which establishes Proposition
G.2 Proof of Proposition
3
Proposition
follows from an application of Weyl's inequality. We have
(A) s or (4)
1 - ( S (o)
since 01 (A) = Mo 2 VM Vnnmax-- Similarly, we have
(A 2 (A) oa 11, 3 (Hoa(E)-VMMApppa 2 MMa(EE),
which establishes the proposition.
G.3 Proof of Proposition 4
We begin by constructing matrices Q and E, via symmetric dilation trick, such that the spectral norms of Q (A A)R and Q EQ coincide, and then apply classical E-nct argument to the spectral norm of Q T EQ, following the proof of Lemma D.1 in [66].
First, we set E = Đ(Ā A), where D is the dilation operator (see Section R11.2 and Q = (Q R). and observe that
Q
A A R (QEQ@Q]
c Q EQu
max ||
where the second equality follows from the Courant-Fischen min-max theorem. Now, let S be an e-net of the d - 1-dimensional unit sphere Sd-1 = {r: ||l|| = 1}. By definition, for any v E Sa-1, there exists some w(v) E S such that || w(v)l22 < E and
Q EQ
max
EQul
||
max
{o
w(v) + w())
Q EQ w(o)+w())
||
S (s² + 28)
Q EQ
+ max WESH
w Q EQw
With = = 1/3, have
Q EQ
ww Q EQu
max
wes!
Now, S can be selected so that its cardinality can be upper bounded by | < 18d (see, for example, Pollard (70).. For a fixed w E S//21 we let z Qu and note that since S0-11 c sa-1 that ||=||2 1, and
h(M+1) ((M+1))(M+1)
n
nd
w Q EQw
Σ
Σ
= 2
>
Ciy 24 2n
i=1
j=1
31

---

Now, over the event that entries egj x M log n. for all ر.J, which occurs which overwhelming
probability by Lemma 3 Hoeffding's 's inequality and union bound over S1/3 gives
2/
P Q EQ 2 >1} R
Σ WEST 1/3
P
Q
EQW
> 9
n
nM Σ
Σ
P
City antt
6
WEST
i=1
2t
2 18 exp
(9cM log
m)²
26
= exp dllo(118) (9cM log n)²
Setting & M log 3/2 n gives
Q (A-A) A) R = =1[ [Q]l 2 M 3/2
log n,
completing the proof.
G.4 Proof of Proposition 5
Denote the singular value decomposition of Ü by REEN, where = == diag((1..... Ea), and let
W == 52, S2 The principal angles [0.]æ1 between the column spaces of Ü and Ü are defined by
& = co,) and by the Wedin sinO theorem, we have
1 | W
I| max & max |1 cos 0:1 S max =E[d]
11 cos- 0:l
2
(E(d]
iE[d]
(21)
A
Mnnmax
n.maa
Vnimmax
maxsin's 0,
S
(A) da+ + (A))2
M
82
We apply the Wedin sinO theorem again to obtain a bound which we will require later:
J00 - 00 00 v /V - vV 1/2
|si θ (Ü, 11,
v sin e V.)ll
|| A
aa (A) 1
ndmmax
5
We start by showing that We now establish a bound on - 8/1-
vTV/
V
- = arg max
(UTO- -UÜ-VAL VVV
x1101
Tj (ŪTÔ-
arg max
x 11 1,=1
+ , (1 (UD-VV))
S arg max
(1 + 5.
x: i=
Σ T ( +8 (0O0 $-1 x
arg max
21 1,ي=1
= (
32

---

and then we employ the decomposition
V ( 8-1
U SV V -- S-1
1) _ i 0/8-1 8-1
U (y-A) 1) v V (^ AA-A)
U ( (A) (v AAA v) S - + -π- VV
+ V (A-A) ( (ū-00 S- + +V (^ x) UUT
Therefore we have
L00-VVV S A-A A11 (1 ( - vV 1/2 + : - oū/²) S
TTT A AA-A) A v S
+
V A (a U
8
AnAmax
1
M log 3/2 n
Mn
+
(A(A)
AA(A)
nAmax
vM log n
82
5
VnAmax
Combining this with 21. we have
P Vnndmax
V w S
-TV-000
-
w/ 2
+
2
G.5 Proof of Proposition 6
We begin by decomposing WS- SW as
w$ Sw
(W- U 0 S $ V - -VV +
= W U ú 8 0000++++5000
V - w) - ( A- 1)
Λ-A) v
(W- (W-U 5 ( V w) +U (Λ-A) a) (VVT -VVT) A
+0T (A A AV
and therefore we have that
|| |W vWS-Sw/ Sw/
< W U Ú
A + w V || l2
A A
VVT VV + U Λ A
AA- v/1
2
ndmmax
VMndmax + M 2/g n
S M 102 n,
which completes the proof.
G.6 Proof of Proposition 7
key tool in proving Proposition is a theorem due [1 providing entrywise eigenvector bounds for random matrices. The original statement is given for the eigenvectors of symmetric random matrices
33

---

with row and column-wis independence. We state generalisation for the singular vectors of rectan- gular matrices with bloek-wise independence structure. The extension block-wise independence structure has been handled in [22 (see Proposition 2. (b) of that paper), although the exposition of the results in this paper more complicated. For this reason, we choose state the result due to [[]] with this generalisation, which can be seen by following through the relevant parts of their proof.
Lemma 8 (A slight generalisation of Theorem 2.1 of (.]. Let Mo be an Th X real-values + T22 and let No == 01 (EM1) /oa ôu = (a(M)M)) Ja+1 (EMo). Suppose there exists some 7 > and function yp R 1R+ which is continuous and non-decreasing on R4, with sp(0) 0 and ppx)//a non-increusing on R,,, such that the following conditions hold:
random matrix. Define no
2 and T2 / f nu. Define
B1 (Incoherence). |MMol2,00 V [EM]
2,00 500
B2 (Block-wise independence). Assume that for any h E [mz]. [n2]. there exists NV C [n1] and N/ C [ma). such that the kth row of Mo independent of the columns {i J & NG and the fth column of Mo is independent of the rows {i : & N2). Let mg maxk,e {INE V [W and assume mo S 6g-
B3 (Spectral norm concentration). Kn max (7, (()) S 1 and IP ( EMoll 12 > > S 1/o for
some 1/o E (0, 1).
B4. [Row and column concentration] There exists some 7) (0, such that for any matrices
Q E Reixd R E R and E [mi].j] E [m2].
( br /Tuobo }}
IP (Mo EMo) R
s dobs
21
11g
and
P -
EMO), a112
S oobboop
bF
}
171
/nobso
n2p
where boc nu Q11,, V T2l 22 and by (TI
[R|
1/2
Let Uo, Uo (respectively Vo. Vo) be the matrices containing the left (respectively, right) singular vectors corresponding to the leading singular values of Mo and EMg Then, with probability at least 1 - Do 2m,, we have
m Toll2.0 V x2ll Vollz,o S (no + (()) (m) Up 12,,2 2,.00 V V V
+ 7 EM, 2..0 V T2 (EM,) 112,0))885
T1 |,, 2 V V Voll2,00 I [No (no + p (1)) {r + (())) (0)1 (T) Uo V Vollz,e0)
7 x EM, 2,00 V n. | 1l2.00) )/80.
The following is adaptation of Lemma D.2 of [66] (see also Lemma 7 of (I] who showed an analogous result for Bernoulli random variables.
Lemma 9. Let Y; ~ Poisson(A:) (X:) independently for all i = I,...,n,, n, and suppose Q is a deter- ministic matrix. The Qi denote the ith row of Q and set Amax max:elm] ); Then for any 0,
n
(2+ 2+aa)nmaas QQ122
(Y+ Σ
x,)a
2de -nnmaax
1 v log
(v | 2, ||Q][
We omit the proof of Lemma 9 which is identical the proof of Lemma D.2 of [66] with the Bernoulli moment generating function with the Poisson moment generating function.
With these tools to hand, we begin by obtaining a bound on | using Lemma with Mg = A.
We set := Vnn.max/8 88 and
n
((x) =
8 log(1x))
34

---

First observe that = + nM x nM and TI x M- and x 1, and that Ko = and ón = VM8. B1 holds since |1A12,,0 I vnd X Vndmax since Amax S by Assumption
Using Assumptions 3 and 4 we have
maxax
5log(6/ /
8 olg log' - 6,
M x
log n
x log
and therefore B2 holds.
B3
holds from Proposition and observing that by Assumption 3
Kg max (7, ((7) S.
To that B4 holds, note that each row and column of M- M A) contains independent Poisson
random variables with means not exceeding ndmax /M. Then for Q Raxd R E BeMxx, setting
a = log n/ max in Lemma implies that
A A
R 2,00
M
A
A
R
M
4,*
2,,00
(nd
M + log n))
|
220000
M
/moal 2
1 V
R
m.mmax
|, 2,00
V /nol|Rlll,ss ||
1 log
R.Il
8 |||12,000 2,00 to
Tio R 2,
OF
I Subocck
Veugoo
Similarly, setting a = M log n/nAmax we have
[ (A-A)
2
v Mndmax
Q112,
(Vudin
Q
I 5obocs0
2,00
1 log
Vno Q
0..2000/ / [Q01=
by
which establishes
B4
Having established
B4 we are ready to apply Lemma
M (No ((1) (T/
V VV2
2,00 (T A
VT2
A
12,00)
180
2,00
have
n nAmaxll /d/n d/n
V 12.00 2.00 = -AASSUS-11 22000 < A
-
Ü
S 2
z
2,00
(2)
(22)
ndmaxct /d/n
2 nM ulgn.
MM
where we used Assumption in the final inequality. Therefore
|1200 S V d T2
ni |
nM plog n,
2,00
and we have
0 (1((())
ndmax (p(1
K
1/2
(Σ)
0
and the first term satisfies
VM (No + ((1) 0.00 VT V 2,0 S VMP(1) (T ArllUl22, V Tel 220
n.maax
M
"I
8
Ammax Vnd logn
35

---

To control the second term, we first observe that
T
2,
M 1/2
A
2,00 M +M-1/211 (M-1/
1,20
M-
| 2,00
|A + M-
max sup ri()) (E(00)
d
2 Mn t VMK + MV ddmax log
S ná uóvdd dd
where the final inequality follows from & S logn I v and Amax S 1. Similarly we obtain
T2l A 12 S u mdmma log using (2 and therefore
A
A
2,00
V n2
A
112.00 S ud ddmax log n.
We then have
MM (TI ||
V
A
AAI12.00
12.00 100 S VM.
Vndmar udv ddmax log 1
v M5
UV uV ndd logn
Combining these bounds, we obtain
2
Wv AidAmax log
2,00
We now apply Lemma 8 with M, = =A(m) We set and ((x) as before and verify Assumptions B1
in the same way. By analogous calculations to the above, we obtain the bound
2
pVnddmaa log
mm)
2,00
8
The final bound is shown in the same way, requiring the additional observation that & fr V p(()) S which follows from Assumption and we obtain
m 10,0~ Uoll2. V Vx2l v. Voll2,00
$ [mo (No + ((11) ( + x(>)) ((1) (m Volla, V V222/ Voll2,0) +
T1 / 2,00 V VTz2 ||A
12.0)))) 180
$ ((1) (T1/|U0l2,.0 V m2 Voll2.00) +n
(T) ||A
2,000 V2 A
180
2,00
uv nnddmax logn
G.7 Proof of Proposition 8
We show
15
using simple application of Wedin's 's inequality:
Â A(m) A
A A
- ūū || = | θ ( $ aa (A) od+1 (A)
2
σd (A)
od+1 (A)
P VMnXmax
v /ndmmxx
25
VM6
8
The proof of (16 requires more delicate argument. We apply Wedin's 's theorem to obtain
(m) A)
(A(m) (m) A) (()))
((m)(( (mm))) - T
sin θ (ū(m), 0)|1, 2
2
A
1d+1 A
ad
(23)
36

---

By Weyl's inequality
od A 2 od (A) +
A -Al, 112 1 ~ od (A).
P
and
aa+1 A S xd+1 (A) -1|ā A
2 0d+1 (A)
and therefore
A
Ja-1 (
ao (A) da+1 (A) = VM 1/2
od ( (2)
0 1/2
d+1 (x))
VM.. (24)
oa
We now focus attention on obtaining bound for || (m) A)ū ||.. Let
Nm fm + (e - 1)m, [M]}
The ijth entry of A (m)
(m) (em)-A) A) =
(om)-A) (A(m)- 1[i
mi== E Nm),
and so (m)) write
A is independent of À (m) and hence A A is independent of Ü(7)) We can then
( (m)) - A) T (m)) = E
(Ame- - Amee-Ame) )
Aie A {(m)
+
REN i=1
=: 11 + C2-
The (square root of the) first term is easily bounded as
1/2
Ä
A
T
2,20
2,00
IA
A
(((m)
2
2,00
HVV HVVdAAAan Addmax log
M Mn.max
5
8
and to bound the second term, we employ Lemma to obtain
1/2 2 M log" n
(m)
+ / /Maxx logn (m)
2,00
F
M log" n
(mm)
+ V M1 log
((())
2,00
2,00
v M Mdmax logn
(m))
2,00
Hv Wwdddmaa log n
M Io
8
2/2 log
3/2
3/2 n
8
where we used Assumption the third inequality. Therefore
(A(m)- - A) T (m))
< >(AAm) - A) T [(m)
1/2 P
2 p log 3/2 n
< < + 6 2
1/2
37

---

Similar analysis yields an analogous bound for (A(m) and 24 we have
A) (m)) I|2, and combining this, with 23
(A(m)
A)
(m))
v v - A) v(mm) ||
(A(m)
( (06m)) -Ū0
Fd
A
td+1 (A)
3/2 log
3/2
2
82
which establishes the proposition.
38


Reviews:
[{"review_id": "lVn4jLRAW3", "reviewer": "Reviewer_SnVH", "summary": "The authors propose an approach for learning time-varying node embeddings from continuous-time dynamic network data, which consist of a set of instantaneous timestamped relational events between nodes (e.g., messages from one social media user to another). Their proposed approach learns a projection that minimizes reconstruction error of the pairwise intensities between nodes and comes with theoretical guarantees on estimation error. They also show that their approach generates embeddings that both preserve network structure at a given time and is temporally coherent. They demonstrate strong empirical performance on simulated data compared to other dynamic network embeddings. Furthermore, they use their approach to analyze a real network data set on face-to-face interactions of primary school students, which is quite enlightening due to the interpretability of their model.\n\n*After rebuttal:* The authors have clarified the one question I had about the meaning of \"inductive\" in their setting. I continue to strongly support the paper.", "questions": "1. The authors mention several times that their approach is inductive, allowing one to obtain a node representation profile outside of the training sample. If the task is to obtain the node representation for the future, how would the Intensity Profile Projection approach handle it? Would it require some data from other nodes at that future time?", "limitations": "Limitations are thoroughly discussed in Section 6. I commend the authors for being very forthcoming with these limitations. I don't view the limitations as weaknesses, because they are mostly limitations that apply to all unsupervised problems.", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "- Proposed approach learns time-varying node embeddings from continuous-time networks with theoretical guarantees, which is among the first, if not the first, in the literature.\n- Proposed embeddings can satisfy two good properties of structure preservation and temporal coherence.\n- Very well written and organized paper that provides highlights of theoretical analysis in the main paper followed by details, including proofs, in the supplementary.", "weaknesses": "- There's a large body of related literature on probabilistic generative models for continuous-time networks using point process models such as Hawkes processes that should be discussed. Many of these models are based on stochastic block models or latent space models and are thus also learning node embeddings. See suggested references below.\n- No quantitative evaluation. This is only a minor weakness in my opinion because I view the main contribution to be theoretical.\n\nTypos and minor issues:\n- Supplementary Section C heading: Visualsation -> Visualisation\n\nReferences:\n- Arastuie, M., Paul, S., & Xu, K. S. (2020). CHIP: A Hawkes process model for continuous-time networks with scalable and consistent estimation. In Advances in Neural Information Processing Systems 33 (pp. 16983-16996).\n- Corneli, M., Latouche, P., & Rossi, F. (2018). Multiple change points detection and clustering in dynamic networks. Statistics and Computing, 28(5), 989-1007. doi:10.1007/s11222-017-9775-1\n- Huang, Z., Soliman, H., Paul, S., & Xu, K. S. (2022). A mutually exciting latent space Hawkes process model for continuous-time networks. In Proceedings of the 38th Conference on Uncertainty in Artificial Intelligence (Vol. 180, pp. 863-873).\n- Junuthula, R. R., Haghdan, M., Xu, K. S., & Devabhaktuni, V. K. (2019). The Block Point Process Model for continuous-time event-based dynamic networks. In Proceedings of the World Wide Web Conference (pp. 829-839).\n- Matias, C., Rebafka, T., & Villers, F. (2018). A semiparametric extension of the stochastic block model for longitudinal networks. Biometrika, 105(3), 665-680. doi:10.1093/biomet/asy016\n- Yang, J., Rao, V., & Neville, J. (2017). Decoupling homophily and reciprocity with latent space network models. In Proceedings of the Conference on Uncertainty in Artificial Intelligence.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> There's a large body of related literature on probabilistic generative models for continuous-time networks using point process models such as Hawkes processes that should be discussed. Many of these models are based on stochastic block models or latent space models and are thus also learning node embeddings. See suggested references below.\n\nThanks for pointing us to these interesting references, which we will point to in the introduction. (See global response.)\n\n> The authors mention several times that their approach is inductive, allowing one to obtain a node representation profile outside of the training sample. If the task is to obtain the node representation for the future, how would the Intensity Profile Projection approach handle it? Would it require some data from other nodes at that future time?\n\nTo obtain a representation of a node, $i$, at a time, $t$, outside the training sample, one just needs the node's intensity profile vector $\\hat \\Lambda_i(t) = (\\hat \\lambda_{i1}(t),\\ldots,\\hat \\lambda_{in})^\\top$, which is then projected onto the subspace spanned by $\\hat{\\mathbf{U}}_d$. The construction of $\\hat\\Lambda_i(t)$ requires data from the interaction events involving node $i$ around time $t$."}, {"review_id": "9An3VEoXpx", "reviewer": "Reviewer_QsMj", "summary": "To represent the continuous dynamic network, authors provide the framework based on the intensity profile. First, the intensity between nodes is estimated, which produces the intensity profile. Low dimension reduction via SVD is applied on the intensity, and then each node embedding is obtained by the low dimensional subspace.\nAuthor also provide various theoretical analysis about the error bound and the bias-various trade-off. Theoretical analysis as well as empirical analysis on the simulated data demonstrates that the proposed method capture structural preserving and temporally coherent properties. Case study on the real data is conducted to explain the outcome of the proposed framework qualitatively", "questions": "- Figure numbers are all wrong. \n- Section 4 is true for any global subspace projection. Also, both properties could be debatable, not necessarily ideal. For instance, when \\Labmda_{i}(s) = \\Lambda_{i}(t), X_{i}(s) = \\alpha * X_{i}(t) could be more ideal, depending on the interactions among the other nodes. \n- It would be great if authors compare the embedding trajectory for more real data, beyond the specific simulated ones. ", "limitations": "Often, the meaning of each dimension from the SVD decomposition is not clear. This interpretability is not necessarily required for the representation, but this should be addressed when presenting the case study.", "rating": 4, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 2, "strengths": "- Simple but powerful method is proposed\n- Based on the mathematical model, theoretical bound is analyzed and explained.\n- IPP can capture the behavior of a bifurcating block model.", "weaknesses": "- The proposed method is not novel enough. SVD decomposition is a very common technique for the reduction of dimensions, and it often suffers from the long-tailed singular values. \n- Comparison is too limited. The analysis has been made only for the simulated data with figures. More experiments as well as some qualitative results would be great to have.\n- SVD decomposition does not prevent producing negative values at the reconstruction.\n- The proposed projection space is very dependent on the fixed dataset. At least, how to leverage the given embeddings for predictions is not straightforward. Given this, the potential application value is not very clear.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> The proposed method is not novel enough. SVD decomposition is a very common technique for the reduction of dimensions,\n\nThis is unreasonable. Dismissing our algorithm as \"not novel\" because it contains an SVD is like dismissing an optimisation algorithm because it uses SGD. The novelty of the algorithm lies in the whole framework which we propose to learn representations, which employs the SVD as the appropriate tool for optimising the $\\hat R^2$ objective function discussed in Section 2 and Lemma 1. In particular, we point out the careful construction of the matrix $\\hat{\\mathbf{\\Sigma}}$ from the data, which is entirely non-trivial.\n\nThe novelty of the method is exemplified in the following points:\n\n- To our knowledge, IPP is the first provably consistent representation learning algorithm for continuous time dynamic networks in the literature. The uniform error bound in Theorem 1 is the first of its kind for data of this kind. We note additionally that our algorithm is non-parametric, and all existing methods are either model-based (e.g. fit a latent position model) and lack statistical estimation theory, or are entirely heuristic.\n- IPP is the first representation learning algorithm for continuous time dynamic networks which satisfies the desirable properties of \"structure preservation\" and \"temporal coherence\" (defined in Section 4), which are necessary to make even the most basic temporal inferences about node behaviours (e.g. \"do nodes $i$ and $j$ behave in the same way at time $t$?\", and \"does node $i$ change its behaviour between times $t$ and $s$?\")\n\n> It [the SVD] often suffers from the long-tailed singular values.\n\nThis is a good point. In much of the data that we have experimented with (e.g. the school children data in Section 5), the eigenvalues $\\hat{\\mathbf{\\Sigma}}$ decay quickly, however the reviewer correctly observes that in some real-world data they might not. In our theory, this corresponds to Assumptions 2 or 3 being violated. \n\nIn the revision, we will add discussion about how this problem can be identified in practice (e.g. by plotting the singular values) and some possible solutions: \n\n- taking a transformation of the intensity profiles (such as the square root) to temper heavy tails [1].\n- employing a robust subspace estimator, such as robust PCA [2].\n\n[1] Ian Gallagher, Andrew Jones, Anna Bertiger, Carey E. Priebe & Patrick Rubin-Delanchy (2023) Spectral Embedding of Weighted Graphs, Journal of the American Statistical Association.\n[2] Candès, E. J., Li, X., Ma, Y., & Wright, J. (2011). Robust principal component analysis?. Journal of the ACM (JACM).\n\n> Comparison is too limited.. / It would be great if authors compare the embedding trajectory for more real data, beyond the specific simulated ones.\n\nWe believe this simple simulation is sufficient to demonstrate the diverse failure mechanisms of each of the rival methods. We have implemented the rival methods on real data, but simply chose not to include them: with no ground truth we are unable to evaluate them objectively. In the included PDF, we include a plot of the omnibus method (Fig. 1a) and aligned spectral embeddings (Fig. 1b) on real data which, at least visually, seem to fail in the same way as they do in simulated data. We will add these to the supplementary material.\n\n> SVD decomposition does not prevent producing negative values at the reconstruction.\n\nFrom the point of view of representation learning, negative values in the intensity reconstructions are immaterial, since they are never used directly. If one desires intensity estimates which are non-negative, negative intensity values can be set to zero, which can only lower the reconstruction error.\n\n>  The proposed projection space is very dependent on the fixed dataset. At least, how to leverage the given embeddings for predictions is not straightforward. Given this, the potential application value is not very clear.\n\nDue to space constraints, we have mainly focused on the application of the learned representations to unsupervised learning (e.g. clustering, trend estimation, etc.) in the real data section and discussion. However, the structure preservation and temporal coherence properties of the embeddings make them particularly suited to dynamic prediction tasks. Our theory not only suggests that the classifier trained on node representations at time $t$ would perform well at predicting unlabelled nodes at that time (see e.g. [3]) but the \"temporal coherence\" property suggests that same classifier would perform well at predicting node labels at some time $s$, in the future. This is something that to our knowledge no other continuous-time dynamic network embedding algorithm can achieve (see the method comparison in Section 4). We will add a comment about the value of our method for predictive applications to the discussion.\n\n[3] Minh Tang. Daniel L. Sussman. Carey E. Priebe. (2013) Universally consistent vertex classification for latent positions graphs. Ann. Statist.\n\n> Section 4 is true for any global subspace projection. Also, both properties could be debatable, not necessarily ideal. For instance, when $\\Lambda_{i}(s) = \\Lambda_{i}(t), X_{i}(s) = \\alpha * X_{i}(t)$ could be more ideal, depending on the interactions among the other nodes.\n\nThe reviewer correctly identifies that these properties hold for any global subspace projection, and this is precisely the part of the algorithm which we are justifying in this section. Justification of the precise global subspace projection we choose is the subject of Lemma 1 in Section 2.\n\nWe think there may be typo here, and that the reviewer meant: $\\Lambda_{i}(s) = \\alpha * \\Lambda_{i}(t)$ implies $X_{i}(s) = \\alpha * X_{i}(t)$. This is a very interesting suggestion and, remarkably, we *do* in fact guarantee that this property holds! This follows by Lemma 2, combined the linearity of the projection. We will add this observation to the revision."}, {"review_id": "uLfJcortPN", "reviewer": "Reviewer_L5U9", "summary": "The paper presents a framework called Intensity Profile Projection (IPP) for continuous-time representation learning in dynamic networks. The authors aim to address the challenge of capturing temporal dynamics and evolving relationships in dynamic networks with both high statistical precision and interpretability. The model leverages the concept of intensity profiles, which encode the temporal changes and interactions between nodes in a network. The model provides a uniform error bound for learned node representations and preserves a novel \"temporal coherence\" property compared to existing baselines. Empirical results on real-world dynamic network datasets demonstrate that IPP outperforms existing methods in various tasks such, highlighting its ability to capture continuous-time representations and uncover temporal patterns in dynamic networks.", "questions": "1. What would be the time complexity of the proposed method?\n2. How would the model performance be affected by different network topology?", "limitations": "The authors have discussed the limitations of their work.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper introduces the Intensity Profile Projection (IPP) framework, which offers a unique and innovative approach to continuous-time representation learning for dynamic networks. It introduces the concept of intensity profiles and effectively utilizes them to capture temporal dynamics. \n\n2. Theoretical analysis towards the model shows that the model can achieve high statistical precision and preserve interpretability in terms of \"\"temporal coherence\".\n\n3. The paper is in general easy to follow.", "weaknesses": "1. Lack of comparison with state-of-the-art methods: Although the paper claims improved performance over existing methods, it does not provide a comprehensive comparison with some existing continuous models such as GraphODEs[1,2,3,4] which combines neuralODE with GNNs to model network evolution over time.\n\n2. Scalability: The scalability of the IPP framework is not extensively discussed. It would be valuable to address the computational requirements and scalability limitations of the proposed approach, especially when dealing with large-scale dynamic networks.\n\n3. The related work section is too short to provide a comprehensive background of the research topic.\n\n\n[1] Huang, Zijie, Yizhou Sun, and Wei Wang. \"Learning continuous system dynamics from irregularly-sampled partial observations.\" Advances in Neural Information Processing Systems 33 (2020): 16177-16187.\n\n[2] Song Wen, Hao Wang, and Dimitris Metaxas. 2022. Social ODE: Multi-agent Trajectory Forecasting with Neural Ordinary Differential Equations. In Computer Vision–ECCV 2022: 17th European Conference.\n\n[3]Zijie Huang, Yizhou Sun, and Wei Wang. Coupled graph ode for learning interacting system dynamics. In\n401 ACM SIGKDD Conference on Knowledge Discovery and Data Mining, page 705–715, 2021.\n\n[4] Zang, Chengxi, and Fei Wang. \"Neural dynamics on complex networks.\" In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pp. 892-902. 2020.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "\n> Lack of comparison with state-of-the-art methods: Although the paper claims improved performance over existing methods, it does not provide a comprehensive comparison with some existing continuous models such as GraphODEs[1,2,3,4] which combines neuralODE with GNNs to model network evolution over time.\n\nWe thank the reviewer for pointing us to GraphODEs. From what we understand, these methods work on multivariate time series on the nodes of a graph which captures the influence of agent $i$ on the dynamics of agent $j$. Our method works on point processes observed between pairs of of nodes. Therefore, these methods are not directly comparable. We are only aware of one existing method (CLPM in our method comparison) which does continuous-time representation learning with the type of data we consider. This is the reason we additionally compare to a collection of discrete time methods.\n\n> Scalability: The scalability of the IPP framework is not extensively discussed. It would be valuable to address the computational requirements and scalability limitations of the proposed approach, especially when dealing with large-scale dynamic networks. / What would be the time complexity of the proposed method?\n\nThis is a great question which we will address in Section 2.2 after discussing the numerical approximation. One of the advantages of the IPP framework is that it is highly scalable. If the network is sparse, the method is scalable to networks with potentially tens of millions of nodes.\n\nSuppose we use a kernel with finite support of width $2h$, then the expected time complexity of the KDE for a single edge at a single point in time is $O(h \\lambda_{\\max})$. Evaluation of $\\hat{\\mathbf{\\Lambda}}(t)$ is $O(n^2 h \\lambda_{\\max})$ and of the $\\hat{\\mathbf{\\Lambda}}$ is $O(Bn^2 h \\lambda_{\\max})$. Consider the special case discussed in Section 3.1 that $\\lambda_{ij} \\asymp \\rho$ for all $i,j$ and $L \\asymp \\rho L_0$. Assuming $\\mu, d$, and $L_0$ are fixed, our theory suggests choosing $h \\asymp (n\\rho)^{-1/3}$, and if we suppose $n \\rho \\asymp \\log^3$, the sparsest regime our theory allows, then evaluation of $\\hat{\\mathbf{\\Lambda}}$ is $O(B n \\log^2 n)$, i.e. log-linear in $n$.\n\nIn practice the top singular vectors of $\\hat{\\mathbf{\\Lambda}}$ can be using the Augmented Implicitly Restarted Lanczos Bidiagonalization algorithm implemented in the `irlba` package in R, or the `irlbpy` in Python. The time complexity of this algorithm has not been studied theoretically although in practice it can be incredibly fast. For example, the package author performs a simulated experiment in which they compute the first 2 singular vectors of a sparse 10M x 10M  matrix with 1M non-zero entries which takes approximately 6 seconds on a computer with two Intel Xeon CPU E5-2650 processors (16 physical CPU cores) running at 2 GHz equipped with 128 GB of ECC DDR3 RAM (see `https://bwlewis.github.io/irlba/comparison.html`).\n\n> The related work section is too short to provide a comprehensive background of the research topic.\n\nSee global response.\n\n> How would the model performance be affected by different network topology?\n\nOur method is able to capture both homophilic and heterophilic network topologies. Fundamentally, the quality of the learned representations depends on how well the intensity profiles are approximated on a d-dimensional subspace."}, {"review_id": "tcs7IDj3Gt", "reviewer": "Reviewer_fL1b", "summary": "The paper proposes a method to learn low-dimensional continuous-time representations of network nodes, based on the collection of interaction events among them. More precisely, the events are in the form of $(i,j,t)$, where $(i,j)$ is the pair of nodes involved in the interaction event, and $t$ is the occurrence time. The proposed method first estimate the intensity function $\\lambda_{i,j}(t)$ of events between each pair of nodes $(i,j)$ at every time instant $t$, then project the intensities of each node at time $t$ onto a learned lower dimensional subspace to obtain a representation. Theoretical results on the recovery error of the representation is provided. Numerical experiments using real data shows the effectiveness of the proposed method.", "questions": ".", "limitations": ".", "rating": 3, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The paper proposes to estimate the representation of nodes using continuous-time events, which seems to be a novel type of data.", "weaknesses": "I find the presentation of the paper generally vague and hand-wavy. See the following.\n\n1. The introduction is way too high-level. The authors should be more specific about the problem setting in this paper, for example, why we care about dynamic models, continuous-time event data, low-dimensional representation of nodes etc.\n\n2. The related work is not specific. The authors should use a sentence to summarize the contribution of the mentioned papers and explain the difference from your work.\n\n3. Lemma 1 is not correct. $\\widehat U_d$ minimizes the residual sum of squares at $B$ chosen time instants, but not the integrated one. \n\n4. In Section 3, notation part, what is the difference between $\\gg$ and $\\gtrsim$? Also is the universal constant multiplicative or additive?\n\n5. It's not clear what `$\\approx$' means in Section 4. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "\n> The introduction is way too high-level. The authors should be more specific about the problem setting in this paper, for example, why we care about dynamic models, continuous-time event data, low-dimensional representation of nodes etc.\n\nSee global response.\n\n> The related work is not specific. The authors should use a sentence to summarize the contribution of the mentioned papers and explain the difference from your work.\n\nSee global response.\n\n> Lemma 1 is not correct. $\\hat U_d$ minimizes the residual sum of squares at $B$ chosen time instants, but not the integrated one.\n\nLemma 1 is correct as stated. Please see the proof in Section D of the supplementary materials.\n\n> In Section 3, notation part, what is the difference between ≫ and ≳? Also is the universal constant multiplicative or additive?\nThanks for pointing this out, we will make this explicit in the revision. The symbols $\\gg$ and $\\gtrsim$ denote the inequalities $>$ and $\\geq$ that hide multiplicative constants. \n\n> It's not clear what '≈' means in Section 4.\n\nThe \"$\\approx$ \" symbol means \"approximately equal to\" and is deliberately left informal in the descriptions of the \"Structure preserving\" and \"Temporally coherent\" properties. In the paragraph before Lemma 2, $\\hat X_i(s) \\approx \\hat X_j(t)$ is formally defined to mean that $X_i(s) = X_j(t)$, so that the lemma is a mathematically rigorous statement."}]

Author rebuttal:
We thank all reviewers for their time and expertise. We summarise some of the positive comments made by reviewers: "a unique and innovative approach to continuous-time representation learning for dynamic networks", "Simple but powerful method", "among the first, if not the first, in the literature.", "Very well written and organized paper", "The paper is in general easy to follow.", "I commend the authors for being very forthcoming with these limitations".

Based on the reviews, we feel we should better convey the concrete possibilities offered by this new algorithm for important applications such as cyber-security [1], combating human-trafficking [2], fraud and corruption [3]. We briefly mentioned these in the conclusion, but in our revision we will discuss these further at the outset.

In addition, we will expand the 'related work' section to cover the broader literature which puts this work in context, for example on estimation of stochastic block model [4-10] and latent position models for continuous time networks [11-16], spectral methods discrete-time dynamic networks [17-19] (in the absence of existing spectral methods for continuous-time networks), probabilistic modelling of network point processes [20-23] (such as Hawkes processes), neural-network based embedding algorithms [24-26] and relevant surveys [27-31].

[1] Kent, A. D. (2015). Cybersecurity Data Sources for Dynamic Network Research. In Dynamic Networks in Cyber-security. Imperial College Press.

[2] Szekely, P., Knoblock, C. A., Slepicka, J., Philpot, A., Singh, A., Yin, C., Kapoor, D., Natarajan, P., Marcu, D., Knight, K. et al. (2015). Building and using a knowledge graph to combat human trafficking. In The Semantic Web-ISWC 2015: 14th International Semantic Web Conference. Springer.

[3] Microsoft Researcg. (2021, December 9). Revealing the Hidden Structure of Corruption. https://www.microsoft.com/en-us/research/group/societal-resilience/articles/revealing-the-hidden-structure-of-corruption/.

[4] Blundell, C., Beck, J., and Heller, K. A. (2012). Modelling reciprocating relationships with hawkes processes. Advances in Neural Information Processing Systems.

[5] DuBois, C., Butts, C., and Smyth, P. (2013). Stochastic blockmodeling of relational event dynamics. In Artificial intelligence and statistics. PMLR.

[6] Corneli, M., Latouche, P., and Rossi, F. (2016). Block modelling in dynamic networks with non-homogeneous Poisson processes and exact ICL. Social Network Analysis and Mining.

[7] Matias, C., Rebafka, T., and Villers, F. (2018). A semiparametric extension of the stochastic block model for longitudinal networks. Biometrika.

[8] Corneli, M., Latouche, P., & Rossi, F. (2018). Multiple change points detection and clustering in dynamic networks. Statistics and Computing.

[9] Junuthula, R., Haghdan, M., Xu, K. S., and Devabhaktuni, V. (2019). The block point process model for continuous-time event-based dynamic networks. In The world wide web conference.

[10] Arastuie, M., Paul, S., and Xu, K. (2020). CHIP: A Hawkes process model for continuous-time networks with scalable and consistent estimation. Advances in Neural Information Processing Systems.

[11] Durante, D. and Dunson, D. B. (2014). Nonparametric Bayes Dynamic Modelling of Relational Data. Biometrika.

[12] Durante, D., & Dunson, D. B. (2016). Locally Adaptive Dynamic Networks. The Annals of Applied Statistics.

[13] Yang, J., Rao, V., & Neville, J. (2017). Decoupling homophily and reciprocity with latent space network models. In Proceedings of the Conference on Uncertainty in Artificial Intelligence.

[14] Rastelli, R. and Corneli, M. (2021). Continuous latent position models for instantaneous interactions. arXiv preprint arXiv:2103.17146.

[15] Huang, Z., Soliman, H., Paul, S., & Xu, K. S. (2022). A mutually exciting latent space Hawkes process model for continuous-time networks. In Proceedings of the 38th Conference on Uncertainty in Artificial Intelligence.

[16] Artico, I. and Wit, E. (2023). Fast inference of latent space dynamics in huge relational event networks. arXiv preprint arXiv:2303.17460.

[17] Liu, F., Choi, D., Xie, L., and Roeder, K. (2018). Global spectral clustering in dynamic networks. Proceedings of the National Academy of Sciences.

[18] Cape, J. (2021). Spectral analysis of networks with latent space dynamics and signs.

[19] Gallagher, I., Jones, A., & Rubin-Delanchy, P. (2021). Spectral embedding for dynamic networks with stability guarantees. Advances in Neural Information Processing Systems.

[20] Butts, C. T. (2008). A relational event framework for social action. Sociological Methodology.

[21] Vu, D., Hunter, D., Smyth, P., and Asuncion, A. (2011). Continuous-time regression models for longitudinal networks. Advances in Neural Information Processing Systems.

[22] Perry, P. O. and Wolfe, P. J. (2013). Point process modelling for directed interaction networks. Journal of the Royal Statistical Society: SERIES B: Statistical Methodology.

[23] Passino, F. S. and Heard, N. A. (2022). Mutually exciting point process graphs for modeling dynamic networks. Journal of Computational and Graphical Statistics.

[24] Nguyen, G. H., Lee, J. B., Rossi, R. A., Ahmed, N. K., Koh, E., and Kim, S. (2018). Continuous-time dynamic network embeddings. In Companion Proceedings of the Web Conference 2018.

[25] Du, L., Wang, Y., Song, G., Lu, Z., and Wang, J. (2018). Dynamic network embedding: An extended approach for skip-gram based network embedding. In IJCAI.

[26] Xu, D., Ruan, C., Korpeoglu, E., Kumar, S., and Achan, K. (2020). Inductive representation learning on temporal graphs. arXiv preprint arXiv:2002.07962.

[27] Spiliopoulou, M. (2011). Evolution in social networks: A survey. Social network data analytics.

[28] Holme, P. and Saram ̈aki, J. (2012). Temporal networks. Physics reports.

... [29-31] omitted due to character limit.

Meta-review: The paper proposes an approach for learning time-varying node embeddings from continuous-time dynamic network data, which consist of a set of instantaneous timestamped relational events between nodes. This is a very interesting and novel paper that has a great number of things going for it. Here is a summary of the positive aspects:

- The paper introduces the Intensity Profile Projection (IPP) framework, offering an innovative approach to continuous-time representation learning for dynamic networks.
- The model utilizes the concept of intensity profiles to capture temporal dynamics.
- Theoretical analyses indicate that the model achieves high statistical precision and preserves "temporal coherence".
- The IPP can capture the behavior of a bifurcating block model.
- Proposed approach learns time-varying node embeddings from continuous-time networks with theoretical guarantees.
- Embeddings can satisfy properties of structure preservation and temporal coherence.
- The paper is very well-written and organized with details provided in supplementary sections.
- Empirical results on simulated data and real-world datasets showcase the effectiveness of the IPP framework.
- The paper proposes to estimate node representations using continuous-time events, a seemingly novel data type.

The authors provided a very strong response to criticism raised in the reports. 
I encourage the authors to try and incorporate components of the response into the camera ready version of the paper, even in the appendinx if the main paper does not have enough space.
Decision: Accept (poster)