# Spectral embedding for dynamic networks with stability guarantees

Ian Gallagher
University of Bristol, UK
<EMAIL_ADDRESS>
&Andrew Jones
University of Bristol, UK
<EMAIL_ADDRESS>
&Patrick Rubin-Delanchy
University of Bristol, UK
<EMAIL_ADDRESS>

###### Abstract

We consider the problem of embedding a dynamic network, to obtain time-
evolving vector representations of each node, which can then be used to
describe changes in behaviour of individual nodes, communities, or the entire
graph. Given this open-ended remit, we argue that two types of stability in
the spatio-temporal positioning of nodes are desirable: to assign the same
position, up to noise, to nodes behaving similarly at a given time (cross-
sectional stability) and a constant position, up to noise, to a single node
behaving similarly across different times (longitudinal stability). Similarity
in behaviour is defined formally using notions of exchangeability under a
dynamic latent position network model. By showing how this model can be recast
as a multilayer random dot product graph, we demonstrate that unfolded
adjacency spectral embedding satisfies both stability conditions. We also show
how two alternative methods, omnibus and independent spectral embedding,
alternately lack one or the other form of stability.

##  1 Introduction

Consider a dynamic network in which nodes come and go, change abruptly or
evolve, alone or in communities, and form connections accordingly. Our goal is
to find a vector representation,
$\hat{\mathbf{Y}}_{i}^{(t)}\in\mathbb{R}^{d}$, for every node $i$ and time
$t$, which could be used for a diversity of downstream analyses, such as
clustering, time series analysis, classification, model selection and more.
This problem is known as dynamic (or evolutionary) network (or graph)
embedding, and a great number of scalable and empirically successful
techniques have been put forward, with recent surveys by [47, 51]. We
consider, among these, a subset about which it is reasonable to try to
establish a certain statistical guarantee.

The novelty of this paper is _not_ to propose a new procedure. Instead, it is
to demonstrate that an existing procedure, unfolded adjacency spectral
embedding (UASE) [17], has two important stability properties. Given a
sequence of symmetric adjacency matrices
$\mathbf{A}^{(1)},\ldots,\mathbf{A}^{(T)}\in\\{0,1\\}^{n\times n}$, where
$\mathbf{A}^{(t)}_{ij}=1$ if nodes $i$ and $j$ form an edge at time $t$, UASE
computes the rank $d$ matrix factorisation of
$\mathbf{A}:=(\mathbf{A}^{(1)}|\cdots|\mathbf{A}^{(T)})$ to obtain
$\mathbf{A}\approx\hat{\mathbf{X}}\hat{\mathbf{Y}}^{\top}$ using the singular
value decomposition (precise details later). The matrix
$\hat{\mathbf{Y}}\in\mathbb{R}^{nT\times d}$ contains, as rows, the desired
representations
$\hat{\mathbf{Y}}_{1}^{(1)},\ldots,\hat{\mathbf{Y}}_{n}^{(1)},\ldots,\hat{\mathbf{Y}}_{1}^{(T)},\ldots,\hat{\mathbf{Y}}_{n}^{(T)}.$

The original motivation for UASE was to analyse multilayer (or multiplex)
graphs under an appropriate extension of the random dot product graph model
[17]. To evaluate UASE and other procedures on the task of dynamic network
embedding, we instead consider the dynamic latent position model

$\displaystyle\mathbf{A}^{(t)}_{ij}\overset{ind}{\sim}\mathrm{Bernoulli}\left(f\left\\{\mathbf{Z}^{(t)}_{i},\mathbf{Z}^{(t)}_{j}\right\\}\right),$
(1)

for $1\leq i<j\leq n$, $t\in[T]$, where
$\mathbf{Z}_{i}^{(t)}\in\mathbb{R}^{k}$ represents the unknown position of
node $i$ at time $t$, and
$f:\mathbb{R}^{k}\times\mathbb{R}^{k}\rightarrow[0,1]$ is a symmetric
function.

This model is well-established [39, 22, 15, 16, 33, 23, 9, 43, 11, 10], with
recent reviews by [20] and [45], and inference usually proceeds on the basis
of a parametric model for $f$ (e.g. logistic in the latent position distance
[39]) and for the dynamics of $\mathbf{Z}_{i}^{(t)}$ (e.g. a Markov process
[39]). The model also includes the dynamic degree-corrected, mixed-membership
and standard stochastic block models as special cases, which were studied in
[48, 52, 14, 29, 50, 49, 31, 4, 32, 19].

To make statistical sense of UASE under this latent position model, we must
somehow connect its output $\hat{\mathbf{Y}}_{i}^{(t)}$ to
$\mathbf{Z}_{i}^{(t)}$. To this end we construct a canonical representative of
$\mathbf{Z}_{i}^{(t)}$, denoted $\mathbf{Y}_{i}^{(t)}$, that UASE can be seen
to estimate. Although we make some regularity assumptions on $f$ and the
$\mathbf{Z}_{i}^{(t)}$, they are not modelled in an explicit, parametric way,
and UASE can clearly be used in practice without having a specific model in
mind. For example, we do not make a Markovian assumption on the evolution of
$\mathbf{Z}_{i}^{(t)}$ and UASE can be used to uncover periodic behaviours.

The key purpose of imposing a dynamic latent position model is to allow us to
put down certain embedding stability requirements. Using this framework, we
can define precisely what we mean by two nodes, $i$ and $j$, behaving
“similarly” at times $s$ and $t$ respectively. In such cases, ideally we would
have $\hat{\mathbf{Y}}_{i}^{(s)}\approx\hat{\mathbf{Y}}_{j}^{(t)}$. Two
special cases are: to assign the same position, up to noise, to nodes behaving
similarly at a given time (cross-sectional stability) and a constant position,
up to noise, to a single node behaving similarly across different times
(longitudinal stability).

To achieve both cross-sectional and longitudinal stability is generally
elusive. We show that two plausible alternatives, omnibus [26] and independent
embedding of each $\mathbf{A}^{(t)}$, alternately exhibit one form of
stability and not the other. More generally, we find existing procedures [6,
27, 40, 8, 54, 7, 28, 5, 4, 32, 37, 19] tend to trade one type of stability
off against the other, e.g. via user-specified cost functions. As a side-note,
it could be observed that omnibus embedding is not being evaluated on a task
for which it was designed, since in the theory of [26] the graphs are
identically distributed. Because the technique is so different from the
others, we still feel it makes an interesting addition.

Our central contribution is to prove that UASE asymptotically provides both
longitudinal and cross-sectional stability: for two nodes, $i$ and $j$,
behaving similarly at times $s$ and $t$ respectively, we have
$\hat{\mathbf{Y}}_{i}^{(s)}\approx\hat{\mathbf{Y}}_{j}^{(t)}$ and, moreover,
$\hat{\mathbf{Y}}_{i}^{(s)}$ and $\hat{\mathbf{Y}}_{j}^{(t)}$ have
asymptotically equal error distribution. We emphasise that these properties
hold without requiring any sort of global stability — the network could vary
wildly over time but certain nodes still stay fixed. In the asymptotic regime
considered, we have $n\rightarrow\infty$, but $T$ fixed so that, for example,
our results are relevant to the case of two large graphs. The alternative
regime where $T\rightarrow\infty$ grows but $n$ is fixed is not easily handled
by existing theory and to provide constant updates to UASE (or omnibus
embedding) in a streaming context presents significant computational
challenges.

The remainder of this article is structured as follows. Section 2 gives a
pedagogical example demonstrating the cross-sectional and longitudinal
stability of UASE in a two-step dynamic stochastic block model, while
highlighting the instability of omnibus and independent spectral embedding. In
Section 3, we prove a central limit theorem for UASE under a dynamic latent
position model, and demonstrate that the distribution satisfies both stability
conditions. In Section 4, we review the stability of other dynamic network
embedding procedures. Section 5 presents an example of UASE applied to a
dynamic network of social interactions in a French primary school, with a
dynamic community detection example given in the main text and a further
classification example provided in the Appendix. Section 6 concludes.

##  2 Motivating example

Suppose it is of interest to uncover dynamic community structure, including
communities changing, merging or splitting. The dynamic stochastic block model
[52, 50] provides a simple explicit model for this, in which two nodes connect
at a certain point in time with probability only dependent on their current
community membership. Suppose that at times 1 and 2, we have the following
inter-community link probability matrices,

$\mathbf{B}^{(1)}=\left(\begin{matrix}0.08&0.02&0.18&0.10\\\
0.02&0.20&0.04&0.10\\\ 0.18&0.04&0.02&0.02\\\
0.10&0.10&0.02&0.06\end{matrix}\right),\quad\mathbf{B}^{(2)}=\left(\begin{matrix}0.16&0.16&0.04&0.10\\\
0.16&0.16&0.04&0.10\\\ 0.04&0.04&0.09&0.02\\\
0.10&0.10&0.02&0.06\end{matrix}\right).$

At time 1 there are four communities present, for example, a node of community
1 connects with a node of community 3 with probability 0.18. At time 2, this
matrix changes so that communities 1 and 2 have merged, community 3 has moved,
whereas community 4 is unchanged. For simplicity, in this example, the
community membership of each node is fixed across time-steps.

We simulate a dynamic network from this model over these two time steps, on
$n=1000$ nodes, equally divided among the four communities and investigate the
results of three embedding techniques: UASE, omnibus, and independent spectral
embedding, displayed in Figure 1. Note that the different techniques produce
embeddings of different dimensions; UASE embeds into $d=4$ dimensions, omnibus
embedding $\tilde{d}=7$, while independent spectral embedding has $d_{1}=4$
and $d_{2}=3$. For visualisation, we show the leading two dimensions for each
embedding. Given the dynamics described above, we contend that the following
properties would be desirable:

  1. 1.

_Cross-sectional stability_ : The embeddings for communities 1 and 2 at time 2
are close.

  2. 2.

_Longitudinal stability_ : The embeddings for community 4 at times 1 and 2 are
close.

Figure 1: First two dimensions of the embeddings for the adjacency matrices
$\mathbf{A}^{(1)}$ and $\mathbf{A}^{(2)}$ using three different techniques:
UASE, omnibus and separate embedding. The points are coloured according to
true community membership, with the black dots showing the fitted community
centroid and the ellipses a fitted 95% level Gaussian contour.

Figure 1 illustrates that UASE has both properties: the blue and orange point
clouds merge at time 2, and the red point cloud has the same location _and_
shape over the two time points. On the other hand, omnibus embedding only has
longitudinal stability (the blue and orange point clouds don’t merge at time
2), whereas independent spectral embedding only has cross-sectional stability
(the red point cloud is at different locations over the two time points).

##  3 Theoretical background and results

We consider a sequence of random graphs distributed according to a dynamic
latent position model, in which the latent position sequences
$(\mathbf{Z}^{(t)}_{i})_{t\in[T]}$ (where we use the shorthand
$[T]=\\{1,\ldots,T\\}$) are independent of each other, and identically
distributed according to a joint distribution $\mathcal{F}$ on
$\mathcal{Z}^{T}$, for a bounded subset $\mathcal{Z}$ of $\mathbb{R}^{k}$. We
emphasise that for a fixed $i$ the latent positions $\mathbf{Z}^{(t)}_{i}$ and
$\mathbf{Z}^{(t^{\prime})}_{i}$ may have different distributions — we are
imposing exchangeability over the nodes (invariance to node labelling), but
not over time.

By extending the function $f$ to be zero outside of its domain of definition,
we may view it as an element of $L^{2}(\mathbb{R}^{k}\times\mathbb{R}^{k})$,
and consequently may define a compact, self-adjoint operator $A$ on
$L^{2}(\mathbb{R}^{k})$ by setting

$\displaystyle Ag(x)=\int_{\mathbb{R}^{k}}f(x,y)g(y)dy$ (2)

for all $g\in L^{2}(\mathbb{R}^{k})$. The operator $A$ has a—possibly
infinite—sequence of non-zero eigenvalues
$\lambda_{1}\geq\lambda_{2}\geq\ldots$ with corresponding orthonormal
eigenfunctions $u_{1},u_{2},\ldots\in L^{2}(\mathbb{R}^{k})$, such that
$Au_{j}=\lambda_{j}u_{j}$ (for further details, see for example [35]). We
shall make the simplifying assumption that $A$ has a _finite_ number of non-
zero eigenvalues, in which case $f$ admits a canonical eigendecomposition

$\displaystyle
f(\mathbf{x},\mathbf{y})=\sum_{i=1}^{D}\lambda_{i}u_{i}(\mathbf{x})u_{i}(\mathbf{y}),$
(3)

assumed to hold everywhere, where after relabelling we may assume that
$|\lambda_{1}|\geq\ldots\geq|\lambda_{D}|$.

Several families of functions satisfy the finite rank assumption $D<\infty$,
such as the multivariate polynomials [35]. Moreover, several existing
statistical models can be written as dynamic latent position network models in
which $D<\infty$, including the dynamic mixed membership, degree-corrected,
and standard stochastic block models. To assume $D<\infty$, more generally, is
tantamount to a claim that, for large $n$, $\mathbf{A}^{(t)}$ has ‘low’
approximate rank: an overwhelming proportion of its eigenvalues are close to
zero. Large matrices with low approximate rank are routinely encountered
across numerous disciplines, and the study [46] provides a hypothesis for this
“puzzling” general observation. However, given a real network, we might reject
the hypothesis that $f$ has low rank (we must then also reject any type of
stochastic block model), for example on the basis of triangle counts [41]. In
such a setting, we anticipate that UASE is still consistent and stable if $d$
is allowed to grow sufficiently slowly with $n$ and there exist asymptotic
results for adjacency spectral embedding, for the single graph case, showing
convergence in Wasserstein distance under assumptions on eigenvalue decay
which can be related to the smoothness of $f$ [24]. However, even if we could
obtain such results for UASE, they would not be as powerful as those we
present here for the finite rank case, where we show uniform consistency with
asymptotic Gaussian error.

######  Theorem 1.

The adjacency matrices $\mathbf{A}^{(1)},\ldots,\mathbf{A}^{(T)}$ are jointly
distributed according to a multilayer random dot product graph model.

The multilayer random dot product graph model (MRDPG, see [17]) is a multi-
graph analogue of the generalised random dot product graph model [36], and is
characterised by the existence of matrices
$\mathbf{\Lambda}^{(t)}\in\mathbb{R}^{d\times d_{t}}$ and random matrices
$\mathbf{X}\in\mathbb{R}^{n\times d}$ and
$\mathbf{Y}^{(t)}\in\mathbb{R}^{n\times d_{t}}$ (generated according to some
joint distribution $\mathcal{G}$) such that

$\displaystyle\mathbf{A}^{(t)}_{ij}\overset{ind}{\sim}\mathrm{Bernoulli}\bigl{(}\mathbf{X}_{i}\mathbf{\Lambda}^{(t)}\mathbf{Y}^{(t)\top}_{j}\bigr{)}$
(4)

for each $t\in[T]$ and $i,j\in[n]$, where $\mathbf{X}_{i}$ and
$\mathbf{Y}^{(t)}_{j}$ denote the $i$th and $j$th rows of $\mathbf{X}$ and
$\mathbf{Y}^{(t)}$ respectively.

We refer the reader to the supplemental material for full details of the proof
of Theorem 1, but note that given knowledge of the function $f$ and the
underlying distribution $\mathcal{F}$ we can construct explicit maps
$\varphi:\mathbb{R}^{Tk}\to\mathbb{R}^{d}$ and
$\varphi_{t}:\mathbb{R}^{k}\to\mathbb{R}^{d_{t}}$, producing row vectors, and
matrices $\mathbf{\Lambda}^{(t)}\in\mathbb{R}^{d\times d_{t}}$ such that

$\displaystyle
f(\mathbf{Z}^{(t)}_{i},\mathbf{Z}^{(t)}_{j})=\mathbf{X}_{i}\mathbf{\Lambda}^{(t)}\mathbf{Y}_{j}^{(t)\top},$
(5)

where
$\mathbf{X}_{i}=\varphi(\mathbf{Z}_{i}),\mathbf{Y}^{(t)}_{j}=\varphi_{t}(\mathbf{Z}_{j}^{(t)})$
and
$\mathbf{Z}_{i}=(\mathbf{Z}^{(1)}_{i}|\cdots|\mathbf{Z}^{(T)}_{i})\in\mathbb{R}^{Tk}$,
for each $t\in[T]$ and $i,j\in[n]$. Consequently, each Gram matrix

$\displaystyle\mathbf{P}^{(t)}=\bigl{(}f(\mathbf{Z}^{(t)}_{i},\mathbf{Z}^{(t)}_{j})\bigr{)}_{i,j\in[n]}$
(6)

admits a factorisation into a product of low-rank matrices, in which the
matrix whose rows are the vectors $\varphi(\mathbf{Z}_{i})$ appears as a
common factor. The dimensions $d_{t}$ and $d$ are precisely the ranks of the
matrices $\mathbf{P}^{(t)}$ and their concatenation
$\mathbf{P}=(\mathbf{P}^{(1)}|\cdots|\mathbf{P}^{(T)})$ respectively. In the
theory that follows we will assume that the dimension $d$ is known and fixed,
whereas in practice we would usually have to estimate it (for example by using
profile likelihood [55]).

We can extend our model to incorporate a range of sparsity regimes by scaling
the function $f$ by a sparsity factor $\rho_{n}$, which we assume is either
constant and equal to $1$, or else tends to zero as $n$ grows, corresponding
to dense and sparse regimes respectively. When this factor is present, the
maps $\varphi$ and $\varphi_{t}$ are scaled by a factor of $\rho_{n}^{1/2}$.

Realising our model as an MRDPG allows us to make precise statements about the
asymptotic behaviour of the point clouds obtained through UASE. In particular,
it is the existence of a common factor in the decomposition of each of the
Gram matrices $\mathbf{P}^{(t)}$ that gives rise to the stability properties
previously demonstrated in the embeddings $\hat{\mathbf{Y}}^{(t)}$, as this
matrix in a sense acts as an anchor for the individual point clouds.

In [17] UASE is used to obtain _two_ multi-graph embeddings; the algorithm
below retains only one — the ‘dynamic’ component.

Algorithm 1 Unfolded adjacency spectral embedding for dynamic networks

1:input symmetric adjacency matrices
$\mathbf{A}^{(1)},\ldots,\mathbf{A}^{(T)}\in\\{0,1\\}^{n\times n}$, embedding
dimension $d$

2:Form the matrix
$\mathbf{A}=(\mathbf{A}^{(1)}|\cdots|\mathbf{A}^{(T)})\in\\{0,1\\}^{n\times
Tn}$ (column concatenation)

3:Compute the truncated singular value decomposition
$\mathbf{U}_{\mathbf{A}}\mathbf{\Sigma}_{\mathbf{A}}\mathbf{V}_{\mathbf{A}}^{\top}$
of $\mathbf{A}$, where $\mathbf{\Sigma}_{\mathbf{A}}$ contains the $d$ largest
singular values of $\mathbf{A}$, and
$\mathbf{U}_{\mathbf{A}},\mathbf{V}_{\mathbf{A}}$ the corresponding left and
right singular vectors

4:Compute the right embedding
$\hat{\mathbf{Y}}=\mathbf{V}_{\mathbf{A}}\mathbf{\Sigma}_{\mathbf{A}}^{1/2}\in\mathbb{R}^{Tn\times
d}$ and create sub-embeddings $\hat{\mathbf{Y}}^{(t)}\in\mathbb{R}^{n\times
d}$ where
$\hat{\mathbf{Y}}=(\hat{\mathbf{Y}}^{(1)};\ldots;\hat{\mathbf{Y}}^{(T)})$ (row
concatenation)

5:return node embeddings for each time period
$\hat{\mathbf{Y}}^{(1)},\ldots,\hat{\mathbf{Y}}^{(T)}$

Replacing $\mathbf{A}$ with $\mathbf{P}$ in this construction yields the
_noise-free_ embeddings $\tilde{\mathbf{Y}}^{(t)}$, whose rows are known to be
linear transformations of the vectors $\varphi_{t}(\mathbf{Z}^{(t)}_{i})$ (see
the supplemental material for further details). A desirable property of UASE
is that since the matrices $\mathbf{A}^{(t)}$ follow an MRDPG model there
exist known asymptotic distributional results for the embeddings
$\hat{\mathbf{Y}}^{(t)}$. In order to ensure the validity of these results, we
will make the assumption that the sparsity factor $\rho_{n}$ satisfies
$\rho_{n}=\omega(n^{-1}\log^{c}(n))$ for some universal constant $c>1$.

In the results to follow, UASE is shown to be consistent and stable in a
uniform sense, i.e. the maximum error of any position estimate goes to zero,
under the assumption that the average network degree grows polylogarithmically
in $n$. To achieve this for less than logarithmic growth would be impossible,
by any algorithm, because it would violate the information-theoretic sparsity
limit for perfect community recovery under the stochastic block model [1]. In
practice this means that the embedding will have high variance when the graphs
have few edges. Several approaches have been proposed to make single graph
embeddings more robust (e.g. to sparsity or heterogeneous degrees), such as
based on the regularised Laplacian [3] or non-backtracking matrices [21], but
it is an open and interesting question how to extend them to the dynamic graph
setting to achieve both cross-sectional and longitudinal stability.

The first of our distributional results states that after applying an
orthogonal transformation—which leaves the structure of the resulting point
cloud intact—the embedded points $\hat{\mathbf{Y}}^{(t)}_{i}$ converge in the
Euclidean norm to the noise-free embedded points
$\tilde{\mathbf{Y}}^{(t)}_{i}$ as the graph size grows:

######  Proposition 2.

There exists a sequence of orthogonal matrices
$\tilde{\mathbf{W}}\in\mathrm{O}(d)$ such that

$\displaystyle\max_{i\in\\{1,\ldots,n\\}}\bigl{\|}\hat{\mathbf{Y}}^{(t)}_{i}\tilde{\mathbf{W}}-\tilde{\mathbf{Y}}^{(t)}_{i}\bigr{\|}=\mathrm{O}\Bigl{(}\tfrac{\log^{1/2}(n)}{\rho_{n}^{1/2}n^{1/2}}\Bigr{)}$
(7)

with high probability for each $t$.

The matrices $\tilde{\mathbf{W}}$ are unidentifiable in practice, but are
defined to be the solution to the one-mode orthogonal Procrustes problem

$\displaystyle\tilde{\mathbf{W}}=\operatorname*{arg\,min}_{\mathbf{Q}\in\mathrm{O}(d)}\|\mathbf{U}_{\mathbf{A}}\mathbf{Q}-\mathbf{U}_{\mathbf{P}}\|_{F}^{2}+\|\mathbf{V}_{\mathbf{A}}\mathbf{Q}-\mathbf{V}_{\mathbf{P}}\|_{F}^{2},$
(8)

where $\lVert\cdot\rVert^{2}_{F}$ denotes the Frobenius norm. We emphasise
that the matrix $\tilde{\mathbf{W}}$ is _the same_ for each embedding, and so
similar behaviour in the noise-free embeddings $\tilde{\mathbf{Y}}^{(t)}$ is
captured by UASE.

Our second result states that after applying a _second_ orthogonal
transformation the above error converges in distribution to a fixed
multivariate Gaussian distribution:

######  Proposition 3.

Let $\zeta=(\zeta_{1}|\cdots|\zeta_{T})\sim\mathcal{F}$, and for
$\mathbf{z}\in\mathcal{Z}$ define

$\displaystyle\mathbf{\Sigma}_{t}(\mathbf{z})=\left\\{\begin{array}[]{cc}\mathbb{E}_{\zeta}\Bigl{[}f(\mathbf{z},\zeta_{t})\bigl{(}1-f(\mathbf{z},\zeta_{t})\bigr{)}\cdot\varphi(\zeta)^{\top}\varphi(\zeta)\Bigr{]}&\mathrm{if~{}}\rho_{n}=1\\\
\mathbb{E}_{\zeta}\Bigl{[}f(\mathbf{z},\zeta_{t})\cdot\varphi(\zeta)^{\top}\varphi(\zeta)\Bigr{]}&\mathrm{if~{}}\rho_{n}\to
0.\end{array}\right..$ (11)

Then there exists a deterministic matrix $\mathbf{R}_{*}\in\mathbb{R}^{d\times
d}$ and a sequence of orthogonal matrices $\mathbf{W}\in\mathrm{O}(d)$ such
that, given $\mathbf{z}\in\mathcal{Z}$, for all $\mathbf{y}\in\mathbb{R}^{d}$
and for any fixed $i\in[n]$ and $t\in[T]$,

$\displaystyle\mathbb{P}\Bigl{(}n^{1/2}(\hat{\mathbf{Y}}^{(t)}_{i}\tilde{\mathbf{W}}-\tilde{\mathbf{Y}}^{(t)}_{i})\mathbf{W}\leq\mathbf{y}~{}|~{}\mathbf{Z}^{(t)}_{i}=\mathbf{z}\Bigr{)}\to\Phi\bigl{(}\mathbf{y},\mathbf{R}_{*}\mathbf{\Sigma}_{t}(\mathbf{z})\mathbf{R}_{*}^{\top}\bigr{)}.$
(12)

As in our previous results, the matrices $\mathbf{W}$ and $\mathbf{R}_{*}$ can
be explicitly constructed given knowledge of the underlying function $f$ and
distribution $\mathcal{F}$ (again, we defer full details to the supplemental
material). Note that as in Proposition 2, both constructed matrices are common
to _all_ embeddings.

We can now demonstrate one of the key advantages that UASE holds over other
embedding methods, namely that _UASE exhibits both cross-sectional and
longitudinal stability_ , in a sense that we shall now define. We say that two
space-time positions $(\mathbf{z},t)$ and $(\mathbf{z}^{\prime},t^{\prime})$
are _exchangeable_ if
$f(\mathbf{z},\zeta_{t})=f(\mathbf{z}^{\prime},\zeta_{t^{\prime}})$ with
probability one, where $\zeta=(\zeta_{1}|\cdots|\zeta_{T})\sim\mathcal{F}$,
and that the positions are _exchangeable up to degree_ if
$f(\mathbf{z},\zeta_{t})=\alpha f(\mathbf{z}^{\prime},\zeta_{t^{\prime}})$ for
some $\alpha>0$. Equivalently, $(\mathbf{z},t)$ and
$(\mathbf{z}^{\prime},t^{\prime})$ are exchangeable if conditional on
$\mathbf{Z}^{(t)}_{i}=\mathbf{z}$ and
$\mathbf{Z}^{(t^{\prime})}_{j}=\mathbf{z}^{\prime}$ the $i$th row of
$\mathbf{P}^{(t)}$ and $j$th row of $\mathbf{P}^{(t^{\prime})}$ are equal with
probability one.

######  Definition 4.

Given a generic method for dynamic network embedding, with output denoted
$(\hat{\mathbf{Z}}^{(t)}_{i})_{i\in[n];t\in[T]}$, define the following
stability properties:

  1. 1.

_Cross-sectional stability:_ Given exchangeable $(\mathbf{z},t)$ and
$(\mathbf{z}^{\prime},t)$, $\hat{\mathbf{Z}}^{(t)}_{i}$ and
$\hat{\mathbf{Z}}^{(t)}_{j}$ are asymptotically equal, with identical error
distribution, conditional on $\mathbf{Z}^{(t)}_{i}=\mathbf{z}$ and
$\mathbf{Z}^{(t)}_{j}=\mathbf{z}^{\prime}$.

  2. 2.

_Longitudinal stability:_ Given exchangeable $(\mathbf{z},t)$ and
$(\mathbf{z},t^{\prime})$, $\hat{\mathbf{Z}}^{(t)}_{i}$ and
$\hat{\mathbf{Z}}^{(t^{\prime})}_{i}$ are asymptotically equal, with identical
error distribution, conditional on
$\mathbf{Z}^{(t)}_{i}=\mathbf{Z}^{(t^{\prime})}_{i}=\mathbf{z}$.

The following result then shows that UASE exhibits both types of stability:

######  Corollary 5.

Conditional on $\mathbf{Z}^{(t)}_{i}=\mathbf{z}$ and
$\mathbf{Z}^{(t^{\prime})}_{j}=\mathbf{z}^{\prime}$, the following properties
hold:

  1. 1.

If $(\mathbf{z},t)$ and $(\mathbf{z}^{\prime},t^{\prime})$ are exchangeable
then $\hat{\mathbf{Y}}^{(t)}_{i}$ and $\hat{\mathbf{Y}}^{(t^{\prime})}_{j}$
are asymptotically equal, with identical error distribution.

  2. 2.

If $(\mathbf{z},t)$ and $(\mathbf{z}^{\prime},t^{\prime})$ are exchangeable up
to degree then $\hat{\mathbf{Y}}^{(t)}_{i}$ and
$\alpha\hat{\mathbf{Y}}^{(t^{\prime})}_{j}$ are asymptotically equal and,
under a sparse regime, their error distributions are equal up to scale,
satisfying
$\mathbf{\Sigma}_{t}(\mathbf{z})=\alpha\mathbf{\Sigma}_{t^{\prime}}(\mathbf{z}^{\prime})$.

We can gain insight into our theoretical results by applying them in the
context of common statistical models. Under the dynamic stochastic block model
of Section 2, a pair $(\mathbf{z},t)$ and $(\mathbf{z}^{\prime},t^{\prime})$
are exchangeable if and only if the corresponding rows of $\mathbf{B}^{(t)}$
and $\mathbf{B}^{(t^{\prime})}$ are identical. Proposition 3 then predicts
that the point cloud obtained via UASE should decompose into a finite number
of _Gaussian_ clusters, as we observe in Figure 1. Moreover, when combined
with Corollary 5 the equality of the fourth rows of $\mathbf{B}^{(1)}$ and
$\mathbf{B}^{(2)}$ implies that we should expect _identical_ clusters (centre
and shape) corresponding to the fourth community at both time points
(indicating longitudinal stability) and similarly the equality of the first
and second rows of $\mathbf{B}^{(2)}$ tells us that the communities should be
indistinguishable at the second time point (indicating cross-sectional
stability), behaviours we observe in Figure 1. Under degree-correction [18,
29], Corollary 5 indicates that we should at least observe cross-sectional and
longitudinal stability along ‘rays’ representing communities, a behaviour that
is exhibited in Figure 2.

##  4 Comparison

In this section we investigate the stability properties of alternatives to
UASE. Table 1 shows the stability of the three embedding algorithms described
in Section 2 and two wider classes of algorithms, described below. For
alternatives to UASE, it will be considered sufficient to establish whether an
embedding is stable when applied to the Gram matrices
$\mathbf{P}^{(1)},\ldots,\mathbf{P}^{(T)}$. A method found to be unstable in
this noise-free condition is not expected to be stable when the matrices are
replaced by their noisy observations
$\mathbf{A}^{(1)},\ldots,\mathbf{A}^{(T)}$.

Table 1: Classes of dynamic network embedding algorithms: a brief description, the algorithm complexity and its cross-sectional/longitudinal stability. In the three latter classes, one may replace $\mathbf{A}^{(t)}$ with other matrix representations of the graph (e.g. the normalised Laplacian). Further details in main text. Algorithm |  Description | Complexity | Stability
---|---|---|---
UASE [17] |  Embed $\mathbf{A}=(\mathbf{A}^{(1)}|\cdots|\mathbf{A}^{(T)})$ | $\mathrm{O}(dTn^{2})$ | Both
Omnibus [26] |  Embed $\tilde{\mathbf{A}}$; $\tilde{\mathbf{A}}_{s,t}=(\mathbf{A}^{(s)}+\mathbf{A}^{(t)})/2$ | $\mathrm{O}(\tilde{d}T^{2}n^{2})$ | Longitudinal
Independent |  Embed $\mathbf{A}^{(t)}$ | $\mathrm{O}(\sum_{t}d_{t}n^{2})$ | Cross-sectional
Separate embedding |  Embed $\bar{\mathbf{A}}^{(t)}=\sum_{k}w_{k}\mathbf{A}^{(t-k)}$ |  | Neither
[40, 8, 4, 32, 37, 19] |  |  |
Joint embedding |  $\operatorname*{arg\,min}_{\hat{\mathbf{Y}}^{(1)},\ldots,\hat{\mathbf{Y}}^{(T)}}\alpha\sum_{t}\text{CS}(\hat{\mathbf{Y}}^{(t)},\mathbf{A}^{(t)})$ |  | Neither
[27, 7, 54, 5, 28] |  $\qquad+(1-\alpha)\sum_{t}\text{CT}(\hat{\mathbf{Y}}^{(t)},\hat{\mathbf{Y}}^{(t+1)})$ |  |

The omnibus method computes the spectral embedding of the matrix
$\tilde{\mathbf{A}}\in\\{0,1\\}^{nT\times nT}$ where the block
$\tilde{\mathbf{A}}_{s,t}\in\mathbb{R}^{n\times n}$ is given by
$(\mathbf{A}^{(k)}+\mathbf{A}^{(\ell)})/2\in\mathbb{R}^{n\times n}$ and we
denote by $\tilde{\mathbf{P}}\in[0,1]^{nT\times nT}$ the noise-free
counterpart of $\tilde{\mathbf{A}}$. If $(\mathbf{z},t)$ and
$(\mathbf{z},t^{\prime})$ are exchangeable, then, conditional on
$\mathbf{Z}^{(t)}_{i}=\mathbf{Z}^{(t^{\prime})}_{i}=\mathbf{z}$, the rows
$\tilde{\mathbf{P}}_{n(t-1)+i}$ and $\tilde{\mathbf{P}}_{n(t^{\prime}-1)+i}$
are equal. However, if $(\mathbf{z},t)$ and $(\mathbf{z}^{\prime},t)$ are
exchangeable, then, in general,
$\tilde{\mathbf{P}}_{n(t-1)+i}\neq\tilde{\mathbf{P}}_{n(t-1)+j}$ when
$\mathbf{Z}^{(t)}_{i}=\mathbf{z}$ and
$\mathbf{Z}^{(t)}_{j}=\mathbf{z}^{\prime}$. One can demonstrate that two rows
of $\tilde{\mathbf{P}}$ are equal if and only if the corresponding nodes’
embeddings are too. Therefore, omnibus embedding provides longitudinal but not
cross-sectional stability.

Independent adjacency spectral embedding computes the spectral embeddings of
the matrices $\mathbf{A}^{(t)}\in\\{0,1\\}^{n\times n}$. If $(\mathbf{z},t)$
and $(\mathbf{z}^{\prime},t^{\prime})$ are exchangeable, then, conditional on
$\mathbf{Z}^{(t)}_{i}=\mathbf{z}$ and
$\mathbf{Z}^{(t^{\prime})}_{j}=\mathbf{z}^{\prime}$, the rows
$\mathbf{P}^{(t)}_{i}$ and $\mathbf{P}^{(t^{\prime})}_{j}$ are equal. If
$t=t^{\prime}$ the embeddings of nodes $i$ and $j$ are equal, however,
embeddings between different graphs are subject to (possibly indefinite)
orthogonal transformations $\mathbf{Q}^{(t)}$ which, if
$\mathbf{P}^{(t)}\neq\mathbf{P}^{(t^{\prime})}$, differ in a non-trivial way
(i.e. _beyond_ simply reflecting the ambiguity of choosing eigenvectors in the
spectral decomposition of $\mathbf{P}^{(t)}$). Therefore, independent
adjacency spectral embedding provides cross-sectional but not longitudinal
stability. The same arguments extend to other independent embeddings.

Separate embedding covers a collection of embedding techniques separately
applied to time-averaged matrices,
$\bar{\mathbf{A}}^{(t)}=\sum_{k}w_{k}\mathbf{A}^{(t-k)}$ where $w_{k}$ are
non-negative weights, and $\mathbf{A}^{(t)}$ may be replaced by another matrix
representation of the graph such as the normalised Laplacian. The weights may
be constant, e.g. $w_{k}=1/t$ for all $k$ [40, 4], exponential forgetting
factors $w_{k}=(1-\lambda)^{k}$ [8, 19], chosen to produce a sliding window
[32], based on a time series model [37], and more. In general, temporal
smoothing results in two nodes behaving identically at time $t$ being embedded
differently if their past or future behaviours differ, whereas the act of
embedding the matrices separately will result in the same issues of alignment
encountered in independent adjacency spectral embedding. Therefore, those
methods can have neither cross-sectional nor longitudinal stability, except in
special cases where one can contrive to have one but not other (e.g.,
$w_{0}=1$ and $w_{k}=0$, reducing to independent embedding). Joint embedding
techniques generally aim to find an embedding to trade-off two costs: a
‘snapshot cost’ (CS) measuring the goodness-of-fit to the observed
$\mathbf{A}^{(t)}$, and a ‘temporal cost’ (CT) penalising change over time.
These are then combined into a single objective function, for example, as [28]
(where we have replaced the normalised Laplacian by the adjacency matrix):

$\operatorname*{arg\,min}_{\breve{\mathbf{A}}^{(1)},\ldots,\breve{\mathbf{A}}^{(T)}}\alpha\sum_{t=1}^{T}\lVert\mathbf{A}^{(t)}-\breve{\mathbf{A}}^{(t)}\rVert^{2}_{F}+(1-\alpha)\sum_{t=1}^{T-1}\lVert\breve{\mathbf{A}}^{(t)}-\breve{\mathbf{A}}^{(t+1)}\rVert^{2}_{F},$
(13)

where $\alpha\in[0,1]$, subject to a low rank constraint on
$\breve{\mathbf{A}}^{(1)},\ldots,\breve{\mathbf{A}}^{(T)}$, where
$\hat{\mathbf{Y}}^{(t)}$ is the spectral embedding of
$\breve{\mathbf{A}}^{(t)}$. It is easy to see that a change affecting only a
fraction of the nodes will result in a change of all node embeddings,
precluding both cross-sectional and longitudinal stability (again apart from
contrived cases, e.g. when $\alpha=1$).

For the complexity calculations, we assume a dense regime in which the
$k$-truncated singular value decomposition of an $m$-by-$n$ matrix is
$\mathrm{O}(kmn)$ [13]. Independent embedding is more efficient than UASE, on
account of $\sum_{t}d_{t}\leq dT$, but UASE is more efficient than omnibus
embedding, because of the linear versus quadratic growth in $T$, while
$\tilde{d}=\operatorname{rank}(\tilde{\mathbf{P}})$ is often larger than $d$.

##  5 Real data

The Lyon primary school data set shows the social interactions at a French
primary school over two days in October 2009 [44]. The school consisted of 10
teachers and 241 students from five school years, each year divided into two
classes. Face-to-face interactions were detected when radio-frequency
identification devices worn by participants (10 teachers and 232 students gave
consent to be included in the experiment) were in close proximity over an
interval of 20 seconds and recorded as a pair of anonymous identifiers
together with a timestamp. The data are available for download from the
Network Repository website111https://networkrepository.com [34].

A time series of networks was created by binning the data into hour-long
windows over the two days, from 08:00 to 18:00 each day. If at least one
interaction was observed between two people in a particular time window, an
edge was created to connect the two nodes in the corresponding network. This
results in a time series of graphs $\mathbf{A}^{(1)},\ldots,\mathbf{A}^{(20)}$
each with $n=242$ nodes. Where a node is not active in a given time window, it
is still included in the graph as an isolated node. This is compatible with
the theory and method, and the node is embedded to the zero vector at that
time point.

Given the unfolded adjacency matrix
$\mathbf{A}=(\mathbf{A}^{(1)}|\cdots|\mathbf{A}^{(20)})$, an estimated
embedding dimension $\hat{d}=10$ was obtained using profile likelihood [55]
and we construct the embeddings
$\hat{\mathbf{Y}}^{(1)},\ldots,\hat{\mathbf{Y}}^{(20)}\in\mathbb{R}^{n\times
10}$, taking approximately five seconds on a 2017 MacBook Pro. Figure 2 shows
the first two dimensions of this embedding to visualise some of the structure
in the data. Similar plots and discussion for both individual spectral
embedding and omnibus embedding are given in the Appendix.

Figure 2: First two dimensions of the embeddings
$\hat{\mathbf{Y}}^{(1)},\ldots,\hat{\mathbf{Y}}^{(20)}$ of the unfolded
adjacency matrix $\mathbf{A}=(\mathbf{A}^{(1)}|\cdots|\mathbf{A}^{(20)})$. The
colours indicate different school years while the marker type distinguish the
two school classes within each year.

From this plot, we observe clustering of students in the same school class.
For time windows corresponding to classroom time, for example, 09:00–10:00 and
15:00–16:00, the embedding forms rays of points in 10-dimensional space, with
each ray broadly corresponding to a single school class. This is to be
expected under a degree-corrected stochastic block model, and the distance
along the ray is a measure of the node’s activity level [25, 30, 38]. However,
not all time windows exhibit this structure, for example, the different
classes mix more during lunchtimes (time windows 12:00–13:00 and 13:00–14:00).

###  5.1 Clustering

Following recommendations regarding community detection under a degree-
corrected stochastic block model [38], we analyse UASE using spherical
coordinates $\mathbf{\Theta}^{(t)}\in[0,2\pi)^{n\times 9}$, for $t\in[T]$.
Since UASE demonstrates cross-sectional and longitudinal stability, we can
combine the embeddings into a single point cloud
$\mathbf{\Theta}=(\mathbf{\Theta}^{(1)\top}|\cdots|\mathbf{\Theta}^{(T)\top})^{\top}\in\mathbb{R}^{nT\times
9}$ where each point represents a student or teacher in a particular time
window. This allows us to detect people returning to a previous behaviour in
the dynamic network. We fit a Gaussian mixture model with varying covariance
matrices to the non-zero points in $\mathbf{\Theta}$ with 20–50 clusters
increasing in increments of 5, with 50 random initialisations, taking
approximately five minutes on a 2017 MacBook Pro. Using the Bayesian
Information Criterion, we select the best fitting model (30 clusters) and
assign the maximum a posteriori Gaussian cluster membership to each student in
each time window.

Figure 3 shows how students in the ten classes move between these clusters
over time. Each class has one or two clusters unique to it, for example, the
majority of students in class 1A spend their classroom time (as opposed to
break time) assigned to cluster 1 or cluster 25. This highlights the
importance of longitudinal stability in UASE, as we are detecting points in
the embedding returning to some part of latent space.

Figure 3: Bar chart showing the Gaussian cluster assignment of school classes
over time. The height of each coloured bar represents the proportion of
students, in that class and at that time, assigned to the corresponding
Gaussian cluster, the total available height representing 100%. If the
coloured bars do not sum to the full available height, the difference
represents the proportion of inactive students. For legibility, only bars
representing over 35% of the class are labelled with the cluster number.

There are also instances of multiple school classes being assigned the same
cluster at the same time period, for example, on the morning of day 1, classes
5A and 5B are mainly in cluster 28 suggesting they are having a joint lesson,
and we see this behaviour again on day 2 with classes 3A and 3B in cluster 20.
In the lunchtime periods, particularly on day 1, the younger students (classes
1A–2B) mingle to form a larger cluster, as do the older students (classes
4A–5B), potentially explained by the cafeteria needing two sittings for lunch
for space reasons [44]. This highlights the importance of cross-sectional
stability in UASE, as it allows the grouping of nodes behaving similarly in a
specific time window, irrespective of their potentially different past and
future behaviours.

##  6 Conclusion

We prove that an existing procedure, UASE, allows dynamic network embedding
with longitudinal and cross-sectional stability guarantees. These properties
make a range of subsequent spatio-temporal analyses possible using ‘off-the-
shelf’ techniques for clustering, time series analysis, classification and
more.

## Funding Transparency Statement

Andrew Jones and Patrick Rubin-Delanchy’s research was supported by an Alan
Turing Institute fellowship. Ian Gallagher’s research was supported by an
EPSRC PhD studentship.

## References

  * Abbe, [2017] Abbe, E. (2017).  Community detection and stochastic block models: recent developments.  The Journal of Machine Learning Research, 18(1):6446–6531.
  * Agterberg et al., [2020] Agterberg, J., Tang, M., and Priebe, C. E. (2020).  On two distinct sources of nonidentifiability in latent position random graph models.  arXiv preprint arXiv:2003.14250.
  * Amini et al., [2013] Amini, A. A., Chen, A., Bickel, P. J., and Levina, E. (2013).  Pseudo-likelihood methods for community detection in large sparse networks.  The Annals of Statistics, 41(4):2097–2122.
  * Bhattacharyya and Chatterjee, [2018] Bhattacharyya, S. and Chatterjee, S. (2018).  Spectral clustering for multiple sparse networks: I.  arXiv preprint arXiv:1805.10594.
  * Chen and Li, [2018] Chen, H. and Li, J. (2018).  Exploiting structural and temporal evolution in dynamic link prediction.  In Proceedings of the 27th ACM International Conference on Information and Knowledge Management, pages 427–436.
  * Chi et al., [2007] Chi, Y., Song, X., Zhou, D., Hino, K., and Tseng, B. L. (2007).  Evolutionary spectral clustering by incorporating temporal smoothness.  In Proceedings of the 13th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 153–162.
  * Deng et al., [2016] Deng, D., Shahabi, C., Demiryurek, U., Zhu, L., Yu, R., and Liu, Y. (2016).  Latent space model for road networks to predict time-varying traffic.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 1525–1534.
  * Dunlavy et al., [2011] Dunlavy, D. M., Kolda, T. G., and Acar, E. (2011).  Temporal link prediction using matrix and tensor factorizations.  ACM Transactions on Knowledge Discovery from Data (TKDD), 5(2):1–27.
  * Durante and Dunson, [2014] Durante, D. and Dunson, D. B. (2014).  Nonparametric Bayes dynamic modelling of relational data.  Biometrika, 101(4):883–898.
  * Durante et al., [2017] Durante, D., Mukherjee, N., and Steorts, R. C. (2017).  Bayesian learning of dynamic multilayer networks.  The Journal of Machine Learning Research, 18(1):1414–1442.
  * Friel et al., [2016] Friel, N., Rastelli, R., Wyse, J., and Raftery, A. E. (2016).  Interlocking directorates in Irish companies using a latent space model for bipartite networks.  Proceedings of the National Academy of Sciences, 113(24):6629–6634.
  * Gao et al., [2011] Gao, S., Denoyer, L., and Gallinari, P. (2011).  Temporal link prediction by integrating content and structure information.  In Proceedings of the 20th ACM international conference on Information and knowledge management, pages 1169–1174.
  * Gu and Eisenstat, [1996] Gu, M. and Eisenstat, S. C. (1996).  Efficient algorithms for computing a strong rank-revealing qr factorization.  SIAM Journal on Scientific Computing, 17(4):848–869.
  * Ho et al., [2011] Ho, Q., Song, L., and Xing, E. (2011).  Evolving cluster mixed-membership blockmodel for time-evolving networks.  In Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics, pages 342–350. JMLR Workshop and Conference Proceedings.
  * Hoff, [2011] Hoff, P. D. (2011).  Hierarchical multilinear models for multiway data.  Computational Statistics & Data Analysis, 55(1):530–543.
  * Hoff et al., [2011] Hoff, P. D. et al. (2011).  Separable covariance arrays via the Tucker product, with applications to multivariate relational data.  Bayesian Analysis, 6(2):179–196.
  * Jones and Rubin-Delanchy, [2021] Jones, A. and Rubin-Delanchy, P. (2021).  The multilayer random dot product graph.  arXiv preprint arXiv:2007.10455.
  * Karrer and Newman, [2011] Karrer, B. and Newman, M. E. J. (2011).  Stochastic blockmodels and community structure in networks.  Phys. Rev. E (3), 83(1):016107, 10.
  * Keriven and Vaiter, [2020] Keriven, N. and Vaiter, S. (2020).  Sparse and smooth: improved guarantees for spectral clustering in the dynamic stochastic block model.  arXiv preprint arXiv:2002.02892.
  * Kim et al., [2018] Kim, B., Lee, K. H., Xue, L., and Niu, X. (2018).  A review of dynamic network models with latent variables.  Statistics surveys, 12:105.
  * Krzakala et al., [2013] Krzakala, F., Moore, C., Mossel, E., Neeman, J., Sly, A., Zdeborová, L., and Zhang, P. (2013).  Spectral redemption in clustering sparse networks.  Proceedings of the National Academy of Sciences, 110(52):20935–20940.
  * Lee and Priebe, [2011] Lee, N. H. and Priebe, C. E. (2011).  A latent process model for time series of attributed random graphs.  Statistical inference for stochastic processes, 14(3):231–253.
  * Lee et al., [2013] Lee, N. H., Yoder, J., Tang, M., and Priebe, C. E. (2013).  On latent position inference from doubly stochastic messaging activities.  Multiscale Modeling & Simulation, 11(3):683–718.
  * Lei, [2021] Lei, J. (2021).  Network representation using graph root distributions.  The Annals of Statistics, 49(2):745–768.
  * Lei and Rinaldo, [2015] Lei, J. and Rinaldo, A. (2015).  Consistency of spectral clustering in stochastic block models.  Ann. Statist., 43(1):215–237.
  * Levin et al., [2017] Levin, K., Athreya, A., Tang, M., Lyzinski, V., Park, Y., and Priebe, C. E. (2017).  A central limit theorem for an omnibus embedding of multiple random graphs and implications for multiscale network inference.  arXiv preprint arXiv:1705.09355.
  * Lin et al., [2008] Lin, Y.-R., Chi, Y., Zhu, S., Sundaram, H., and Tseng, B. L. (2008).  Facetnet: a framework for analyzing communities and their evolutions in dynamic networks.  In Proceedings of the 17th international conference on World Wide Web, pages 685–694.
  * Liu et al., [2018] Liu, F., Choi, D., Xie, L., and Roeder, K. (2018).  Global spectral clustering in dynamic networks.  Proceedings of the National Academy of Sciences, 115(5):927–932.
  * Liu et al., [2014] Liu, S., Wang, S., and Krishnan, R. (2014).  Persistent community detection in dynamic social networks.  In Pacific-Asia Conference on Knowledge Discovery and Data Mining, pages 78–89. Springer.
  * Lyzinski et al., [2014] Lyzinski, V., Sussman, D. L., Tang, M., Athreya, A., and Priebe, C. E. (2014).  Perfect clustering for stochastic blockmodel graphs via adjacency spectral embedding.  Electron. J. Stat., 8(2):2905–2922.
  * Matias and Miele, [2015] Matias, C. and Miele, V. (2015).  Statistical clustering of temporal networks through a dynamic stochastic block model.  arXiv preprint arXiv:1506.07464.
  * Pensky et al., [2019] Pensky, M., Zhang, T., et al. (2019).  Spectral clustering in the dynamic stochastic block model.  Electronic Journal of Statistics, 13(1):678–709.
  * Robinson and Priebe, [2012] Robinson, L. F. and Priebe, C. E. (2012).  Detecting time-dependent structure in network data via a new class of latent process models.  arXiv preprint arXiv:1212.3587.
  * Rossi and Ahmed, [2015] Rossi, R. A. and Ahmed, N. K. (2015).  The network data repository with interactive graph analytics and visualization.  In AAAI.
  * Rubin-Delanchy, [2020] Rubin-Delanchy, P. (2020).  Manifold structure in graph embeddings.  In Proceedings of the Thirty-fourth Conference on Neural Information Processing Systems.
  * Rubin-Delanchy et al., [2020] Rubin-Delanchy, P., Cape, J., Priebe, C. E., and Tang, M. (2020).  A statistical interpretation of spectral embedding: the generalised random dot product graph.  arXiv preprint arXiv:1709.05506v3.
  * Sanna Passino et al., [2019] Sanna Passino, F., Bertiger, A. S., Neil, J. C., and Heard, N. A. (2019).  Link prediction in dynamic networks using random dot product graphs.  arXiv preprint arXiv:1912.10419.
  * Sanna Passino et al., [2020] Sanna Passino, F., Heard, N. A., and Rubin-Delanchy, P. (2020).  Spectral clustering on spherical coordinates under the degree-corrected stochastic blockmodel.  arXiv e-prints, pages arXiv–2011.
  * Sarkar and Moore, [2005] Sarkar, P. and Moore, A. W. (2005).  Dynamic social network analysis using latent space models.  ACM SigKDD Explorations Newsletter, 7(2):31–40.
  * Scheinerman and Tucker, [2010] Scheinerman, E. R. and Tucker, K. (2010).  Modeling graphs using dot product representations.  Computational statistics, 25(1):1–16.
  * Seshadhri et al., [2020] Seshadhri, C., Sharma, A., Stolman, A., and Goel, A. (2020).  The impossibility of low-rank representations for triangle-rich complex networks.  Proceedings of the National Academy of Sciences, 117(11):5631–5637.
  * Seto et al., [2015] Seto, S., Zhang, W., and Zhou, Y. (2015).  Multivariate time series classification using dynamic time warping template selection for human activity recognition.  In 2015 IEEE Symposium Series on Computational Intelligence, pages 1399–1406. IEEE.
  * Sewell and Chen, [2015] Sewell, D. K. and Chen, Y. (2015).  Latent space models for dynamic networks.  Journal of the American Statistical Association, 110(512):1646–1657.
  * Stehlé et al., [2011] Stehlé, J., Voirin, N., Barrat, A., Cattuto, C., Isella, L., Pinton, J.-F., Quaggiotto, M., Van den Broeck, W., Régis, C., Lina, B., et al. (2011).  High-resolution measurements of face-to-face contact patterns in a primary school.  PloS one, 6(8):e23176.
  * Turnbull, [2020] Turnbull, K. (2020).  Advancements in latent space network modelling.  PhD thesis, Lancaster University.
  * Udell and Townsend, [2019] Udell, M. and Townsend, A. (2019).  Why are big data matrices approximately low rank?  SIAM Journal on Mathematics of Data Science, 1(1):144–160.
  * Xie et al., [2020] Xie, Y., Li, C., Yu, B., Zhang, C., and Tang, Z. (2020).  A survey on dynamic network embedding.  arXiv preprint arXiv:2006.08093.
  * Xing et al., [2010] Xing, E. P., Fu, W., Song, L., et al. (2010).  A state-space mixed membership blockmodel for dynamic network tomography.  Annals of Applied Statistics, 4(2):535–566.
  * Xu, [2015] Xu, K. (2015).  Stochastic block transition models for dynamic networks.  In Artificial Intelligence and Statistics, pages 1079–1087. PMLR.
  * Xu and Hero, [2014] Xu, K. S. and Hero, A. O. (2014).  Dynamic stochastic blockmodels for time-evolving social networks.  IEEE Journal of Selected Topics in Signal Processing, 8(4):552–562.
  * Xue et al., [2021] Xue, G., Zhong, M., Li, J., Chen, J., Zhai, C., and Kong, R. (2021).  Dynamic network embedding survey.  arXiv preprint arXiv:2103.15447.
  * Yang et al., [2011] Yang, T., Chi, Y., Zhu, S., Gong, Y., and Jin, R. (2011).  Detecting communities and their evolutions in dynamic social networks—a Bayesian approach.  Machine learning, 82(2):157–189.
  * Zhao et al., [2017] Zhao, B., Lu, H., Chen, S., Liu, J., and Wu, D. (2017).  Convolutional neural networks for time series classification.  Journal of Systems Engineering and Electronics, 28(1):162–169.
  * Zhu et al., [2016] Zhu, L., Guo, D., Yin, J., Ver Steeg, G., and Galstyan, A. (2016).  Scalable temporal latent space inference for link prediction in dynamic social networks.  IEEE Transactions on Knowledge and Data Engineering, 28(10):2765–2777.
  * Zhu and Ghodsi, [2006] Zhu, M. and Ghodsi, A. (2006).  Automatic dimensionality selection from the scree plot via the use of profile likelihood.  Computational Statistics & Data Analysis, 51(2):918–930.

##  Appendix A Appendix

Python notebooks producing the figures of this paper are available at
https://github.com/iggallagher/Dynamic-Network-Embedding.

### Proof of Theorem 1

Given the factorisation

$\displaystyle
f(\mathbf{x},\mathbf{y})=\sum_{i=1}^{D}\lambda_{i}u_{i}(\mathbf{x})u_{i}(\mathbf{y}),$
(14)

define a map $\phi:\mathcal{Z}\to\mathbb{R}^{D}$ by setting the $i$th
coordinate of $\phi(\mathbf{z})$ to be $|\lambda_{i}|^{1/2}u_{i}(\mathbf{z})$,
so that for any $\mathbf{x},\mathbf{y}\in\mathcal{Z}$ we have
$f(\mathbf{x},\mathbf{y})=\phi(\mathbf{x})^{\top}\mathbf{I}_{p,q}\phi(\mathbf{y})$
(where $\mathbf{I}_{p,q}$ is the diagonal matrix whose entries are the signs
of the eigenvalues $\lambda_{i}$). Let $\mathcal{F}^{*}$ be the joint
distribution on $\mathbb{R}^{TD}$ obtained by first assigning a random vector
$\zeta=(\zeta_{1}|\cdots|\zeta_{T})$ via $\mathcal{F}$ and then applying the
map $\phi$ to each of the components $\zeta_{t}$, and let
$\mathcal{F}^{*}_{1},\ldots,\mathcal{F}^{*}_{T}$ denote the corresponding
marginal distributions on $\mathbb{R}^{D}$. Given $\xi\sim\mathcal{F}^{*}$ and
$\xi_{t}\sim\mathcal{F}^{*}_{t}$, define the second moment matrices
$\Delta=\mathbb{E}[\xi\xi^{\top}]\in\mathbb{R}^{TD\times TD}$ and
$\Delta_{t}=\mathbb{E}[\xi_{t}\xi_{t}^{\top}]\in\mathbb{R}^{D\times D}$, and
let $r=\mathrm{rank}(\Delta)$ and $r_{t}=\mathrm{rank}(\Delta_{t})$.

Let $\mathbf{M}\in\mathbb{R}^{r\times TD}$ be a matrix whose rows form a basis
of $\mathrm{supp}(\mathcal{F}^{*})$, and similarly let
$\mathbf{N}_{t}\in\mathbb{R}^{r_{t}\times D}$ be a matrix whose rows form a
basis of $\mathrm{supp}(\mathcal{F}^{*}_{t})$. Let
$\mathbf{N}=\mathrm{diag}(\mathbf{N}_{1},\ldots,\mathbf{N}_{T})$, and define
the matrix
$\mathbf{\Pi}:=\mathbf{M}\mathbf{D}\mathbf{N}^{\top}\in\mathbb{R}^{r\times(r_{1}+\ldots+r_{T})}$,
where $\mathbf{D}=\mathrm{diag}(\mathbf{I}_{p,q},\ldots,\mathbf{I}_{p,q})$.
Construct the singular value decomposition
$\mathbf{\Pi}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$, with
$\mathbf{U}\in\mathrm{O}(r\times d)$, $\mathbf{\Sigma}\in\mathbb{R}^{d\times
d}$ and $\mathbf{V}\in\mathrm{O}((r_{1}+\cdots+r_{T})\times d)$, and
$d=\mathrm{rank}(\mathbf{\Pi})$. Writing
$\mathbf{V}=(\mathbf{V}_{1}|\cdots|\mathbf{V}_{T})$, where
$\mathbf{V}_{t}\in\mathbb{R}^{r_{t}\times d}$ has rank $d_{r}$, we can then
construct the singular value decompositions
$\mathbf{V}_{t}=\mathbf{U}_{t}\mathbf{\Sigma}_{t}\mathbf{W}_{t}^{\top}$, with
$\mathbf{U}_{t}\in\mathrm{O}(r_{t}\times d_{t})$,
$\mathbf{\Sigma}_{t}\in\mathbb{R}^{d_{t}\times d_{t}}$ and
$\mathbf{W}_{t}\in\mathrm{O}(d\times d_{t})$, where
$d_{t}=\mathrm{rank}(\mathbf{V}_{t})$.

Define
$\mathbf{\Lambda}=\mathbf{\Sigma}\mathbf{W}\in\mathbb{R}^{d\times(d_{1}+\cdots+d_{T})}$,
where
$\mathbf{W}=(\mathbf{W}_{1}\mathbf{\Sigma}_{1}|\cdots|\mathbf{W}_{T}\mathbf{\Sigma}_{T})$,
and let $L:\mathbb{R}^{TD}\to\mathbb{R}^{d}$ be the linear map sending the
$i$th row of $\mathbf{M}$ to the $i$th row of $\mathbf{U}$, and similarly let
$L_{t}:\mathbb{R}^{D}\to\mathbb{R}^{d_{t}}$ be the linear map sending the
$i$th row of $\mathbf{N}_{t}$ to the $i$th row of $\mathbf{U}_{t}$. Finally,
let $\varphi:\mathbb{R}^{Tk}\to\mathbb{R}^{d}$ and
$\varphi_{t}:\mathbb{R}^{k}\to\mathbb{R}^{d_{t}}$ be the maps satisfying
$\varphi(\mathbf{z})=L\bigl{(}(\phi(\mathbf{z}^{(1)})|\cdots|\phi(\mathbf{z}^{(T)})\bigr{)}$
and
$\varphi_{t}(\mathbf{z}^{(t)})=L_{t}\bigl{(}\phi(\mathbf{z}^{(t)})\bigr{)}$
for any
$\mathbf{z}=(\mathbf{z}^{(1)}|\cdots|\mathbf{z}^{(T)})\in\mathbb{R}^{Tk}$.

Then, setting $\mathcal{G}$ to be the joint distribution on
$\mathbb{R}^{d}\times\mathbb{R}^{d_{1}}\times\cdots\mathbb{R}^{d_{T}}$
obtained by first assigning a random vector
$\zeta=(\zeta_{1}|\cdots|\zeta_{T})$ via $\mathcal{F}$ and then sending this
to the tuple
$\bigr{(}\varphi(\zeta),\varphi_{1}(\zeta_{1}),\ldots,\varphi_{t}(\zeta_{t})\bigr{)}$
and letting $\mathbf{X}_{i}=\varphi\bigl{(}\mathbf{Z}_{i}\bigr{)}$ and
$\mathbf{Y}_{i}^{(t)}=\varphi_{t}\bigl{(}\mathbf{Z}^{(t)}_{i}\bigr{)}$, we
find that
$(\mathbf{A},\mathbf{X},\mathbf{Y})\sim\mathrm{MRDPG}(\mathcal{G},\mathbf{\Lambda})$.
$\square$

### Proof of Proposition 2

This follows directly from Theorem 1 in [17], which states that there exist
sequences of matrices
$\mathbf{R}_{t}=\mathbf{R}_{t}(n)\in\mathbb{R}^{d_{t}\times d}$ and
$\tilde{\mathbf{W}}\in\mathrm{O}(d)$ such that

$\displaystyle\|\hat{\mathbf{Y}}^{(t)}-\mathbf{Y}^{(t)}\mathbf{R}_{t}\tilde{\mathbf{W}}^{\top}\|_{2\to\infty}=\mathrm{O}\Bigl{(}\tfrac{\log^{1/2}(n)}{\rho^{1/2}n^{1/2}}\Bigr{)}$
(15)

almost surely, where the matrices $\mathbf{R}_{t}$ satisfy
$\tilde{\mathbf{Y}}^{(t)}=\mathbf{Y}^{(t)}\mathbf{R}_{\mathrm{Y},t}$, and
$\tilde{\mathbf{W}}\in\mathrm{O}(d)$ is the solution to the one-mode
orthogonal Procrustes problem

$\displaystyle\tilde{\mathbf{W}}=\operatorname*{arg\,min}_{\mathbf{Q}\in\mathrm{O}(d)}\|\mathbf{U}_{\mathbf{A}}\mathbf{Q}-\mathbf{U}_{\mathbf{P}}\|_{F}^{2}+\|\mathbf{V}_{\mathbf{A}}\mathbf{Q}-\mathbf{V}_{\mathbf{P}}\|_{F}^{2},$
(16)

where $\mathbf{A}$ and $\mathbf{P}$ admit the singular value decompositions
$\mathbf{A}=\mathbf{U}_{\mathbf{A}}\mathbf{\Sigma}_{\mathbf{A}}\mathbf{V}_{\mathbf{A}}^{\top}+\mathbf{U}_{\mathbf{A},\perp}\mathbf{\Sigma}_{\mathbf{A},\perp}\mathbf{V}_{\mathbf{A},\perp}^{\top}$
and
$\mathbf{P}=\mathbf{U}_{\mathbf{P}}\mathbf{\Sigma}_{\mathbf{P}}\mathbf{V}_{\mathbf{P}}^{\top}$
respectively. Observing that the $2$-to-infinity norm of a matrix is known to
be equivalent to its maximum Euclidean row norm (and consequently is invariant
under orthogonal transformations) gives the desired result. $\square$

### Proof of Proposition 3

From the proof of Theorem 1 in [17], we find that

$\displaystyle
n^{1/2}(\hat{\mathbf{Y}}^{(t)}\tilde{\mathbf{W}}-\tilde{\mathbf{Y}}^{(t)})=n^{1/2}(\mathbf{A}^{(t)}-\mathbf{P}^{(t)})\mathbf{U}_{\mathbf{P}}\mathbf{\Sigma}_{\mathbf{P}}^{-1/2}+n^{1/2}\mathbf{E},$
(17)

where the residual term $\mathbf{E}$ satisfies
$\|n^{1/2}\mathbf{E}\|_{2\to\infty}\to 0$. We can rewrite this as

$\displaystyle
n^{1/2}(\hat{\mathbf{Y}}^{(t)}\tilde{\mathbf{W}}-\tilde{\mathbf{Y}}^{(t)})=n^{1/2}(\mathbf{A}^{(t)}-\mathbf{P}^{(t)})\mathbf{X}\mathbf{L}\mathbf{\Sigma}_{\mathbf{P}}^{-1}+n^{1/2}\mathbf{E},$
(18)

where $\mathbf{L}\in\mathrm{GL}(d)$ (the general linear group of invertible
$d\times d$ matrices) satisfies $\mathbf{X}=\tilde{\mathbf{X}}\mathbf{L}$,
which is known to exist by Proposition 16 of [17].

We begin by showing that there exists a sequence of orthogonal matrices
$\mathbf{W}\in\mathrm{O}(d)$ and a fixed matrix
$\tilde{\mathbf{L}}\in\mathbb{R}^{d\times d}$ such that
$\mathbf{L}\mathbf{W}_{*}\to\tilde{\mathbf{L}}$, for which we adapt the
arguments of Theorem 1 and Corollary 2 of [2].

To begin with, note that the mapping
$v\mapsto\mathbf{X}(\mathbf{X}^{\top}\mathbf{X})^{-1/2}v$ sends the
eigenvectors of the matrix
$(\mathbf{X}^{\top}\mathbf{X})^{1/2}\mathbf{\Lambda}\mathbf{Y}^{\top}\mathbf{Y}\mathbf{\Lambda}^{\top}(\mathbf{X}^{\top}\mathbf{X})^{1/2}$
to the left singular vectors of $\mathbf{P}$, since if $v$ is such an
eigenvector (with corresponding eigenvalue $\lambda$) then

$\displaystyle\mathbf{P}\mathbf{P}^{\top}\mathbf{X}(\mathbf{X}^{\top}\mathbf{X})^{-1/2}v$
$\displaystyle=\mathbf{X}\mathbf{\Lambda}\mathbf{Y}^{\top}\mathbf{Y}\mathbf{\Lambda}^{\top}(\mathbf{X}^{\top}\mathbf{X})^{1/2}v$
(19)
$\displaystyle=\mathbf{X}(\mathbf{X}^{\top}\mathbf{X})^{-1/2}(\mathbf{X}^{\top}\mathbf{X})^{1/2}\mathbf{\Lambda}\mathbf{Y}^{\top}\mathbf{Y}\mathbf{\Lambda}^{\top}(\mathbf{X}^{\top}\mathbf{X})^{1/2}v$
(20) $\displaystyle=\lambda\mathbf{X}(\mathbf{X}^{\top}\mathbf{X})^{-1/2}v$
(21)

as required. Thus we may write
$\mathbf{U}_{\mathbf{P}}=\mathbf{X}(\mathbf{X}^{\top}\mathbf{X})^{-1/2}\mathbf{V}$,
where $\mathbf{V}$ is a matrix of eigenvectors of
$(\mathbf{X}^{\top}\mathbf{X})^{1/2}\mathbf{\Lambda}\mathbf{Y}^{\top}\mathbf{Y}\mathbf{\Lambda}^{\top}(\mathbf{X}^{\top}\mathbf{X})^{1/2}$,
and consequently observe that

$\displaystyle\mathbf{L}$
$\displaystyle=(\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{X}_{\mathbf{P}}$
(22)
$\displaystyle=(\mathbf{X}^{\top}\mathbf{X})^{-1}\mathbf{X}^{\top}\mathbf{X}(\mathbf{X}^{\top}\mathbf{X})^{-1/2}\mathbf{V}\mathbf{\Sigma}_{\mathbf{P}}^{1/2}$
(23)
$\displaystyle=(\mathbf{X}^{\top}\mathbf{X})^{-1/2}\mathbf{V}\mathbf{\Sigma}_{\mathbf{P}}^{1/2}$
(24)
$\displaystyle=\Bigl{(}\tfrac{\mathbf{X}^{\top}\mathbf{X}}{n}\Bigr{)}^{-1/2}\mathbf{V}\Bigl{(}\tfrac{\mathbf{\Sigma}_{\mathbf{P}}}{n}\Bigr{)}^{1/2}.$
(25)

Let $\mathcal{G}_{\mathrm{X}}$ and $\mathcal{G}_{\mathrm{Y},t}$ denote the
marginal distributions of $\mathcal{G}$, and define the second moment matrices
$\Delta_{\mathrm{X}}=\mathbb{E}[\xi\xi^{\top}]$ and
$\Delta_{\mathrm{Y},t}=\mathbb{E}[\xi_{t}\xi_{t}^{\top}]$, where
$\xi\sim\mathcal{G}_{\mathrm{X}}$ and $\xi_{t}\sim\mathcal{G}_{\mathrm{Y},t}$,
and let
$\Delta_{\mathrm{Y}}=\mathrm{diag}(\Delta_{\mathrm{Y},1},\ldots,\Delta_{\mathrm{Y},T})$.
Then the law of large numbers tells us that the first and last terms in (25)
converge to $(\rho_{n}\Delta_{\mathrm{X}})^{-1/2}$ and
$(\rho_{n}\tilde{\mathbf{\Sigma}})^{1/2}$ respectively, where
$\tilde{\mathbf{\Sigma}}$ is the diagonal matrix whose entries are the square
roots of the eigenvalues of
$\Delta_{\mathrm{X}}^{1/2}\mathbf{\Lambda}\Delta_{\mathrm{Y}}\mathbf{\Lambda}^{\top}\Delta_{\mathrm{X}}^{1/2}$
(see, for example, Proposition 7 of [17]). Note that $\mathbf{V}$ is also the
matrix of eigenvectors of

$\displaystyle\Bigl{(}\tfrac{\mathbf{X}^{\top}\mathbf{X}}{n}\Bigr{)}^{1/2}\Bigl{(}\tfrac{\mathbf{\Lambda}\mathbf{Y}^{\top}\mathbf{Y}\mathbf{\Lambda}^{\top}}{n}\Bigr{)}\Bigl{(}\tfrac{\mathbf{X}^{\top}\mathbf{X}}{n}\Bigr{)}^{1/2}$
(26)

which converges to
$\Delta_{\mathrm{X}}^{1/2}\mathbf{\Lambda}\Delta_{\mathrm{Y}}\mathbf{\Lambda}^{\top}\Delta_{\mathrm{X}}^{1/2}$
by the law of large numbers. Consequently, for each distinct eigenvalue of
$\Delta_{\mathrm{X}}^{1/2}\mathbf{\Lambda}\Delta_{\mathrm{Y}}\mathbf{\Lambda}^{\top}\Delta_{\mathrm{X}}^{1/2}$
we may apply the Davis-Kahan theorem to find that the principal angles between
the resulting eigenspace and the subspace spanned by the corresponding columns
of $\mathbf{V}$ vanish, and thus $\mathbf{V}$ converges to
$\tilde{\mathbf{V}}$ up to some block-orthogonal transformation
$\mathbf{W}\in\mathrm{O}(d)$, where $\tilde{\mathbf{V}}$ is a fixed matrix of
eigenvectors of
$\Delta_{\mathrm{X}}^{1/2}\mathbf{\Lambda}\Delta_{\mathrm{Y}}\mathbf{\Lambda}^{\top}\Delta_{\mathrm{X}}^{1/2}$.
Since $\mathbf{W}$ by definition commutes with $\tilde{\mathbf{\Sigma}}$, we
find that

$\displaystyle\mathbf{L}\mathbf{W}\to\Delta_{\mathrm{X}}^{-1/2}\tilde{\mathbf{V}}\tilde{\mathbf{\Sigma}}^{1/2}$
(27)

as required.

Multiplying (18) by $\mathbf{W}$, we find that

$\displaystyle
n^{1/2}(\hat{\mathbf{Y}}^{(t)}_{i}\tilde{\mathbf{W}}-\tilde{\mathbf{Y}}^{(t)}_{i})\mathbf{W}\approx
n\rho_{n}\Biggl{[}\frac{1}{n^{1/2}\rho_{n}}(\mathbf{A}^{(t)}-\mathbf{P}^{(t)})\mathbf{X}\Biggr{]}_{i}\mathbf{L}\mathbf{\Sigma}_{\mathbf{P}}^{-1}\mathbf{W}$
(28)

and note that the term
$n\rho_{n}\mathbf{L}\mathbf{\Sigma}_{\mathbf{P}}^{-1}\mathbf{W}$ converges to
$\tilde{\mathbf{L}}\tilde{\mathbf{\Sigma}}^{-1}$ from our previous discussion.
Moreover,

$\displaystyle\Biggl{[}\frac{1}{n^{1/2}\rho_{n}}(\mathbf{A}^{(t)}-\mathbf{P}^{(t)})\mathbf{X}\Biggr{]}_{i}$
$\displaystyle=\frac{1}{n^{1/2}\rho_{n}}\sum_{j=1}^{n}(\mathbf{A}^{(t)}_{ij}-\mathbf{P}^{(t)}_{ij})\mathbf{X}_{j}$
(29)
$\displaystyle=\frac{1}{(n\rho_{n})^{1/2}}\sum_{j=1}^{n}(\mathbf{A}^{(t)}_{ij}-\mathbf{P}^{(t)}_{ij})\varphi(\mathbf{Z}_{j}).$
(30)

Conditional on $\mathbf{Z}_{i}^{(t)}=\mathbf{z}$, we have
$\mathbf{P}^{(t)}_{ij}=\rho_{n}f(\mathbf{z},\mathbf{Z}^{(t)}_{j})$, and so the
sum in (30) is a scaled sum of $n-1$ independent, identically distributed
zero-mean random variables, each with covariance matrix

$\displaystyle\mathbb{E}_{\zeta}\Bigl{[}f(\mathbf{z},\zeta_{t})\bigl{(}1-\rho_{n}f(\mathbf{z},\zeta_{t})\bigr{)}\cdot\varphi(\zeta)^{\top}\varphi(\zeta)\Bigr{]}$
(31)

where $\zeta\sim\mathcal{F}$, from which the result follows by setting
$\mathbf{R}_{*}=\tilde{\mathbf{\Sigma}}^{-1}\tilde{\mathbf{L}}^{\top}$ and
applying the multivariate versions of the central limit theorem and Slutsky’s
theorem. $\square$

### Proof of Corollary 5

Note that for any $t\in[T]$ the equality
$\tilde{\mathbf{X}}\tilde{\mathbf{Y}}^{(t)\top}=\mathbf{P}^{(t)}$ holds, and
so
$\tilde{\mathbf{Y}}^{(t)}=\mathbf{P}^{(t)}\tilde{\mathbf{X}}(\tilde{\mathbf{X}}^{\top}\tilde{\mathbf{X}})^{-1}$
(since $\mathbf{P}^{(t)}$ is symmetric). Consequently, for any $i\in[n]$ we
find that
$\tilde{\mathbf{Y}}^{(t)}_{i}=(\tilde{\mathbf{X}}^{\top}\tilde{\mathbf{X}})^{-1}\tilde{\mathbf{X}}^{\top}\mathbf{P}^{(t)}_{i}$,
where $\mathbf{P}^{(t)}_{i}$ denotes the $i$th row of $\mathbf{P}^{(t)}$
(again due to symmetry of $\mathbf{P}^{(t)}$). Thus for any exchangeable pair
$(\mathbf{z},t)$ and $(\mathbf{z}^{\prime},t^{\prime})$, if
$\mathbf{Z}^{(t)}_{i}=\mathbf{z}$ and
$\mathbf{Z}^{(t^{\prime})}_{j}=\mathbf{z}^{\prime}$ then the equality of the
rows $\mathbf{P}^{(t)}_{i}$ and $\mathbf{P}^{(t^{\prime})}_{j}$ implies the
equality of $\tilde{\mathbf{Y}}^{(t)}_{i}$ and
$\tilde{\mathbf{Y}}^{(t^{\prime})}_{j}$. Moreover, from the definition of
exchangeability it is clear that the matrices
$\mathbf{\Sigma}_{t}(\mathbf{z})$ and
$\mathbf{\Sigma}_{t^{\prime}}(\mathbf{z}^{\prime})$ present in the limiting
distributions for $\hat{\mathbf{Y}}^{(t)}_{i}$ and
$\hat{\mathbf{Y}}^{(t^{\prime})}_{j}$ in Proposition 3 are equal, and since
the matrices $\mathbf{W}$, $\tilde{\mathbf{W}}$ and $\mathbf{R}_{*}$ are
independent of $t$ and $t^{\prime}$ equality of the full covariance matrices
follows.

An analogous argument holds in the case that $(\mathbf{z},t)$ and
$(\mathbf{z}^{\prime},t^{\prime})$ are exchangeable up to degree. $\square$

### Lyon primary school data: other embeddings

In this section, we show the spectral embeddings for the Lyon primary school
data for two other embedding algorithms: independent adjacency spectral
embedding and omnibus embedding.

First, for the adjacency matrices $\mathbf{A}^{(1)},\ldots,\mathbf{A}^{(20)}$,
we construct the adjacency spectral embeddings
$\hat{\mathbf{Y}}^{(1)},\ldots,\hat{\mathbf{Y}}^{(20)}\in\mathbb{R}^{n\times
10}$, using the same dimension as the equivalent UASE for this data. This
takes approximately 0.2 seconds on a 2017 MacBook Pro. Figure 4 shows the
first two dimensions of this embedding to visualise some of the structure in
the data.

Figure 4: First two dimensions of the embeddings
$\hat{\mathbf{Y}}^{(1)},\ldots,\hat{\mathbf{Y}}^{(20)}$ of the adjacency
matrices $\mathbf{A}^{(1)},\ldots,\mathbf{A}^{(20)}$. The colours indicate
different school years while the marker type distinguishes the two school
classes within each year.

From this plot, we see that there is no longitudinal stability using
individual adjacency spectral embedding. The classes shown in the leading two
dimensions are not the same between time periods. Therefore, we cannot make
any inference about how the behaviour of students changes over time.

Secondly, we construct the omnibus matrix
$\tilde{\mathbf{A}}\in\mathbb{R}^{nT\times nT}$, where the $n$-by-$n$ block of
the matrix corresponding to times $s,t$ is given by
$\tilde{\mathbf{A}}_{s,t}=(\mathbf{A}^{(s)}+\mathbf{A}^{(t)})/2\in\mathbb{R}^{n\times
n}$. We construct the spectral embedding of $\tilde{\mathbf{A}}$ into
$\tilde{d}=10$ dimensions (as with the other embeddings), to obtain
$\hat{\mathbf{Y}}\in\mathbb{R}^{nT\times 10}$. This is divided to get the
omnibus embedding for each time period,
$\hat{\mathbf{Y}}=(\hat{\mathbf{Y}}^{(1)}|\cdots|\hat{\mathbf{Y}}^{(20)})$.
This takes approximately 30 seconds on a 2017 MacBook Pro. Figure 5 shows the
first two dimensions of this embedding to visualise some of the structure in
the data.

Figure 5: First two dimensions of the embeddings
$\hat{\mathbf{Y}}^{(1)},\ldots,\hat{\mathbf{Y}}^{(20)}$ of the omnibus matrix
$\tilde{\mathbf{A}}$. The colours indicate different school years while the
marker type distinguish the two school classes within each year.

The rough temporal alignment of the point clouds shows the benefits of omnibus
embedding in providing longitudinal stability. However, the lunchtime periods
(days 1 and 2, 12:00-14:00) best demonstrate the effect of a lack of cross-
sectional stability. Some of the teachers are embedded among what we believe
are largely spurious clusters of students during these lunch periods.
Empirically, teachers interact much less frequently with students in these
periods. Moreover, such a sharp classroom-wise clustering at lunchtime is not
present in the independent embedding and seems inconsistent with personal
experience. We believe these effects are due to the averaging of past and
future behaviours inherent in omnibus embedding, which causes cross-sectional
instability.

### Lyon primary school data: classification

An alternative use of UASE is to analyse the trajectory of each node through
time, in embedded space. Because of the longitudinal stability of UASE,
standard multivariate time series analysis techniques can be used to detect
trends and seasonal behaviour. These, in turn, could enable latent position
forecasting and, from this, link prediction. As background, in such analyses,
the time series model is usually incorporated in the embedding process, for
example, via a Markov model for the communities [5, 12, 54], or a seasonal
autoregressive integrated process on the adjacency matrices [37].

In this section, we instead consider the task of time series classification.
Given the trajectory for each student, the goal is to predict their school
class. Given a arbitrary time series, this is often done using dynamic time
warping [42] or convolutional neural networks [53] to allow for misalignment
in time. However, in this simple example with fixed classroom times, we simply
fit a random forest classifier with 100 trees to the concatenation of the
spectral embeddings
$(\mathbf{\Theta}^{(1)}|\cdots|\mathbf{\Theta}^{(T)})\in\mathbb{R}^{n\times
9T}$, each using five randomly selected features. The 10-fold cross-validation
accuracy is $0.983\pm 0.035$, taking approximately two seconds on a 2017
MacBook Pro.

When classifying a time series in this way, auto-correlation makes feature
importance harder to measure. Nevertheless, certain features appear _not_ to
be important for classification, in particular, we find that the two lunchtime
periods (time windows 12:00-13:00 and 13:00-14:00) are not useful, confirming
what is shown by those spectral embeddings in Figure 2.

Using this classifier, we can predict the class of the 10 teachers. While we
do not expect them to behave exactly as the students, we might hope to match a
teacher to their class. Figure 6 shows the proportion of random forest trees
classifying each teacher to each school class. Without truth data it is
impossible to know if these labels are correct, but the most likely
classification assigns exactly one teacher to each class.

Figure 6: Heat map showing the proportion of random forest trees assigning the
spectral embedding trajectory for each teacher to the ten school classes
1A–5B. The colour represents the school class, matching the colours used in
Figure 2, where larger proportions are more opaque. Missing values represent a
proportion of $0$.