Rawya Zreik

Pierre Latouche

Charles Bouveyron

The Dynamic Random Subgraph Model for the Clustering of Evolving Networks

Keywords: State space model, variational inference, variational expectation maximization, maritime data

In recent years, many clustering methods have been proposed to extract information from networks. The principle is to look for groups of vertices with homogenous connection profiles. Most of these techniques are suitable for static networks, that is to say, not taking into account the temporal dimension. This work is motivated by the need of analyzing evolving networks where a decomposition of the networks into subgraphs is given. Therefore, in this paper, we consider the random subgraph model (RSM) which was proposed recently to model networks through latent clusters built within known partitions. Using a state space model to characterize the cluster proportions, RSM is then extended in order to deal with dynamic networks. We call the latter the dynamic random subgraph model (dRSM). A variational expectation maximization (VEM) algorithm is proposed to perform inference. We show that the variational approximations lead to an update step which involves a new state space model from which the parameters along with the hidden states can be estimated using the standard Kalman filter and Rauch-Tung-Striebel (RTS) smoother. Simulated data sets are considered to assess the proposed methodology. Finally, dRSM along with the corresponding VEM algorithm are applied to an original maritime network built from printed Lloyd's voyage records.

Introduction

Network analysis has become a mature discipline, since the original work of Moreno (1934), which is no longer limited to sociology and is now applied in many areas such as biology (Albert and Barabási, 2002;Barabási and Oltvai, 2004;Palla et al., 2005), geography (Ducruet, 2013) or history (Rossi et al., 2014). The growing interest in network analysis is explained partly by the strong presence of this type of data in the digital world, and by recent advances in the modeling and the processing of these data. The clustering methods allow in particular clusters of vertices sharing homogeneous connection profiles to be uncovered. Most methods look for specific structures, so called communities, which exhibit a transitivity property such that nodes of the same community are more likely to be connected (Hofman and Wiggins, 2008). A popular approach for community discovering, though asymptotically biased (Bickel and Chen, 2009), is based on the modularity score given by Girvan and Newman (2002). Alternative methods usually rely on the latent position cluster model (LPCM) of Handcock et al. (2007) which assumes that the links between the vertices depend on their positions in a social latent space.

The stochastic block model (SBM) (Wang and Wong, 1987;Nowicki and Snijders, 2001) is a flexible random graph model which can also characterize communities, but not only. It is based on a probabilistic generalization of the method applied by White et al. (1976) on Sampson's famous monastery (Fienberg and Wasserman, 1981). The SBM model assumes that each vertex belongs to a latent group, and that the probability of connection between a pair of vertices depends exclusively on their group. Because no assumption is made on the connection probabilities, various types of structures of vertices can be taken into account. While SBM was originally developed to analyze mainly binary networks, many extensions have been proposed since to deal for instance with valued edges (Mariadassou et al., 2010) or to take into account prior information (Zanghi et al., 2010;Matias and Robin, 2014). In particular, the random subgraph model (RSM) of Jernite et al. (2014) aims at modeling categorical edges using prior knowledge of a partition of the network into subgraphs. These known subgraphs are assumed to be made of latent clusters which have to be inferred. The vertices are then connected with a probability depending only on the subgraphs whereas the edge type is assumed to be sampled conditionally on the latent groups. This model was applied in the original paper to analyze a historical network in merovingian Gaul. Note that other extensions of SBM have focused on looking for overlapping clusters (Airoldi et al., 2008;Latouche et al., 2011). The inference of SBM like models is usually done using variational expectation maximization (VEM) (Daudin et al., 2008), variational Bayes EM (VBEM) (Latouche et al., 2012), or Gibbs sampling (Nowicki and Snijders, 2001). Moreover, we emphasize that various strategies have been derived to estimates the number of corresponding clusters using model selection criteria (Daudin et al., 2008;Latouche et al., 2012), allocation sampler (Mc Daid et al., 2013), greedy search (Côme and Latouche, 2015), or non parametric schemes (Kemp et al., 2006).

Recently, a few attempts have been made to extend the models mentioned previously in order to deal with dynamic networks. The main idea consists in introducing temporal processes in order to characterize the temporal evolution of nodes and edges through time. Thus, Yang et al. (2011) proposed a dynamic version of SBM allowing a node to switch its class at time t + 1 depending on its state at time t. The switching probabilities are all characterized by a transition matrix. The alternative extension for SBM of Xu and Hero III (2013) focuses on modeling the temporal changes through a state space model and relies on the Kalman filter for inference. Contrary to Yang et al. (2011), Xu and Hero III (2013) treated the edge probabilities as time varying parameters.

In parallel, the mixed membership SBM (MMSBM) of Airoldi et al. (2008), capable of characterizing overlapping clusters, was adapted to deal with dynamic networks by Xing et al. (2010), Ho et al. (2011) and Kim and Leskovec (2013). Moreover, Sarkar and Moore (2005) derived a dynamic version of the LPCM model of Handcock et al. (2007) keeping the transitivity property that nodes which are close in a social latent space should be more likely to connect. Finally, we would like to highlight the work of Dubois et al. (2013) and Heaukulani and Ghahramani (2013). In Dubois et al. (2013) a non homogeneous Poisson process is considered. Thus, contrary to most clustering models for dynamic networks, a continuous time period is taken into account and events, i.e. the creation or removal of an edge, occur one at a time. While models usually focus on modeling the dynamic of networks through the evolution of their latent structures, Heaukulani and Ghahramani (2013) extended the dynamic latent feature model of Foulds et al. (2011) to define how observed social interactions can affect future unobserved latent structures. In the same vein, a dynamic model inspired by SBM was proposed recently by Xu (2015).

In this paper, we aim at modeling dynamic networks with binary or more generally typed edges, for which a partition of the nodes is given. As an example, we will consider an original network, built from printed Lloyd's voyage records and describing maritime flows between ports where the geographical positions of the ports play an important role. The partition was obtained by associating each port to a region according to its geographical position. Figure 1 presents the evolution of network navigations, for 23 years between October 1985 and October 2008. A (given) partition of the nodes is seen here as a decomposition of the network into known subgraphs that we propose to model using unobserved clusters that have to be inferred from the data in practice. Thus, considering a slightly different version of the original RSM model of Jernite et al. (2014) and relying on a state space model as in Xing et al. (2010), we propose a new random graph model for evolving networks that we call the dynamic RSM (dRSM) model. The model focuses on describing the network dynamic by characterizing the evolution of the cluster proportions within the known subgraphs. A logistic transformation is used to link the hidden states and the clusters proportions, as in Blei and Lafferty (2007a); Ahmed and Xing (2007). The inference of the model is done using a VEM algorithm.

The article is organized as follows. In Section 2, we introduce the dRSM model along with an inference procedure in Section 3. Variational techniques are considered and a model selection criterion is derived. Finally, the methodology is tested on simulated data in Section 4 and on the maritime network built from Lloyd's data in Section 5.

The dynamic random subgraph model

This section presents the context of the work and introduces the dRSM model along with the modeling of its dynamic. The joint distribution associated with the model is also detailed. Data extracted from Lloyd's list. The known subgraphs correspond to geographical regions (continents) indicated using colors.

Context and notations

We consider a set of T networks {G (t) } T t=1 , where G (t) is a directed graph observed at time t. Each G (t) is represented by its N × N adjacency matrix X (t) where N denotes the number of nodes. The edge X (t) i,j , describing the relationship between nodes i and j, is assumed to take its values in {0, . . . C} such that X (t) ij = c means that nodes i and j are linked by a relationship of type c at time t and X (t) ij = 0 indicates the absence of relationship between the two nodes at time t. Note that no self loops are considered, i.e. the connection of a node to itself, thus

X (t) ii = 0, ∀ i, t.
Moreover, a partition P of the network into S classes of vertices is assumed to be given. We emphasize that the observed partition induces a decomposition of the graph into subgraphs where each class of vertices corresponds to a specific subgraph. To describe the subgraph membership of each vertex, the variable s is introduced. The variable takes its values in {1, . . . S} and is such that s i indicates the subgraph of vertex i. In some cases, and in order to clarify the equations, we will also consider the indicator variables y is such that y is = 1 if node i is in subgraph s, 0 otherwise. Finally, because the vertex i can only belong to a single subgraph, we have S s=1 y is = 1. Our goal is to cluster at each time t the N nodes into K latent groups with homogeneous connection profiles, i.e. find an estimate of the set Z of latent variables Z (t) ik such that Z (t) ik = 1 if at time t, the node i belongs to the class k, and 0 otherwise. Please note that N , C, P, S and K are all assumed to be constant over time.

The model at each time t

As in the original RSM model, the (known) subgraphs are assumed to be built from K unobserved clusters of vertices, with varying proportions. Thus, each subgraph s has its own mixing proportion vector α

(t) s = (α (t) s1 , . . . , α (t) sK ) where α (t)
sk is the proportion of cluster k in subgraph s at time t and K k=1 α (t) sk = 1, ∀ s, t. The network is then assumed to be generated at each time t as follows.

Each vertex i is first associated to a latent cluster k with a probability depending on its subgraph s i . In practice, the variable

Z (t) i is drawn from a multinomial distribution of parameter α (t) si : Z (t) i ∼ M(1, α (t) si ),
and therefore

K k=1 Z (t) ik = 1. Note that Z (t)
ik = 1 indicates that vertex i belongs to cluster k at time t, 0 otherwise.

On the other hand, the type of link between nodes i and j is assumed to be sampled from a multinomial distribution depending on the latent vectors

Z (t) i and Z (t) j : X (t) ij |Z (t) ik Z (t) jl = 1 ∼ M(1, Π kl ), with Π kl ∈ [0, 1] C+1 and C c=0 Π c kl = 1, ∀k, l. 1 2 4 5 3 Π •• Π •• 8 7 6 9 Π •• Π •• Π •• Figure 2: A dRSM network observed at time t.
The network is made of 9 nodes belonging to S = 2 subgraphs (denoted through the form of the nodes) and split into K = 3 clusters (indicated by the colors). According to the dRSM model, the directed edges between the nodes can be of different types (C = 2 types are considered here). Given the clusters, the presence of an edge depends on the connection probabilities between clusters (Π).

As in the RSM model, and more generally in SBM like models, all vectors Z (t) i are sampled independently, and, conditionally on these membership vectors, the edges are assumed to be independent. Thus, contrary to the original RSM model, the edges depend directly on the latent clusters exclusively, and there is no direct dependency on the subgraphs (see Figure 3). Each edge between a pair (i, j) of vertices does depend on the subgraphs s i and s j , but only through the fact that the edge depends on the latent clusters of the vertices, which themselves depend on the subgraphs. The dependency is indirect while in the original RSM model, the latent clusters along with the subgraphs are all involved in the creation of edges and have different roles. Indeed, the presence or absence of an edge between (i, j) is first drawn from a Bernoulli distribution depending on s i and s j . If an edge is present, the edge type is then sampled depending on the latent clusters. The separation of roles between the latent clusters and the subgraphs was originally motivated by assumptions regarding the nature of the networks analyzed. We do not make such assumptions in this paper. The latent clusters explain both the creation of an edge and its type. Figure 2 presents an example of a dRSM network, observed at time t, made of 9 nodes belonging to 2 subgraphs (denoted through the form of nodes) and split into 3 clusters (indicated by the colors).

Modeling the evolution of random subgraphs

In order to model the evolution of the cluster proportions within the subgraphs through time, a state space model is considered as in Xing et al. (2010). Thus, the latent variable γ (t) s is introduced and a logistic transformation f (•) is used to link the mixing vector α

(t) s with γ (t) s : α (t) s = f (γ (t) s ), such that α (t) sk = f k (γ (t) s ) = exp(γ (t) sk -C(γ (t) s )), ∀s, k, t,
where γ (t) sK = 0 and C(γ

(t) s ) = log( K k=1 exp(γ (t)
sk )). The choice to fix the last component of the vector γ (t) s arbitrarily to 0 is widely used in the literature (see for instance Blei and Lafferty, 2007a;Lafferty and Blei, 2006;Blei and Lafferty, 2007b;Xing et al., 2010) and is due to the bijectivity constraint of this logistic transformation which requires γ (t) s to live in a (K -1) dimensional space since α (t) s has (K -1) degrees of freedom. This induces that γ

(t) sk = log(α (t) sk /α (t)
sK ), ∀s, k, t. In addition, the (K -1) first components of the vector γ (t) s are assumed to be distributed according to a Gaussian distribution with mean Bν (t) and covariance matrix Σ:

γ (t) s\K ∼ N (Bν (t) , Σ), (1)
where γ

(t)

s\K is the vector γ (t)
s without his last component. Both Σ and B are matrices of size (K -1) × (K -1) while ν (t) is a (K -1) dimensional vector. Let us notice that even though the γ (t) s have the same mean in the state-space, they are actually independent and thus play different roles.

The rest of the model now involves a classic state space model for linear dynamic systems. It is defined as follows:

ν (t) = Aν (t-1) + ω ν (1) = µ 0 + u.
The noise terms ω and u are supposed to be Gaussian and independent:

ω ∼ N (0, Φ) u ∼ N (0, V 0 ).
Again, A, Φ and V 0 are matrices of size (K -1) × (K -1) while µ 0 is a (K -1) dimensional vector.

Notice that the state space model for linear dynamic systems may suffer from model identifiability issues and constraints have to be introduced (see for instance Harvey, 1989). In the following, we derive the inference procedure in a general context since different constraints can be considered. In practice, in all the experiments that we carried out, we fixed A, B, and V 0 to be equal to the identity matrix I K-1 and all components of µ 0 to zero.

The model described here has three sets of latent variables (ν = (ν (t) ) t , γ = (γ

(t) s ) st , Z = (Z (t) ik ) ikt )
and is parameterized by θ = (µ 0 , A, B, Φ, V 0 , Σ, Π). Note that all parameters in θ depend neither on time nor subgraphs. This model is called the dynamic random subgraph model (dRSM) in the rest of the document. Figure 3 gives the graphical model for dRSM and Table 1 summarizes the notations used in the model.

At this point, it is possible to see some links and differences between dRSM and dM3SBM (Ho et al., 2011), which is the closest model in the litterature. On the one hand, dRSM and dM3SBM share a common way to model the latent clusters and the temporal dynamic through a state space model. On the other hand, dRSM is able to handle categorical edges, which is a useful feature when working on real-world networks, whereas dM3SBM cannot. In addition, dRSM requires the knowledge of the subgraphs whereas dM3SBM proposes to estimate them. Furthermore, dM3SBM allows the nodes to belong to different clusters. However, allowing to estimate the subgraphs and multi-group belongings may conduce dM3SBM to be a too flexible model and thus to fail in recovering the network structure. Indeed, providing the subgraphs to dRSM allows it to avoid looking for obvious structures such that it can focus on the search of hidden patterns. The comparisons presented in Section 4 seem to confirm this thesis.

X (t) ij Π Z (t) i Z (t) j P γ (t) B, Σ ν (t) µ 0 , A, Φ, ν 0

Joint distribution of dRSM

The dRSM model proposed above is defined by the joint distribution:

p(X, Z, γ, ν|θ) = p(X|Z, Π)p(Z|γ)p(γ \K |B, ν, Σ)p(ν|µ 0 , A, Φ, V 0 ), (2)
where γ \K = (γ

(t) s\K ) st . Moreover p(X|Z, Π) = T t=1 K k,l C c=0 (Π c kl ) N i =j δ(X (t) ij =c)Z (t) ik Z (t) jl , and
p(Z|γ) = T t=1 N i=1 K k=1 f k (γ (t) si ) Z (t) ik = T t=1 K k=1 S s=1 f k (γ (t) s ) N i=1 yisZ (t)
ik .

(

(t) sk = f k (γ (t)
s ) is the proportion of cluster k in the subgraph s at t

Table 1: Summary of the notations used in the paper.

Note that

p(γ \K |B, ν, Σ) = T t=1 S s=1 N (γ (t) s\K ; Bν (t) , Σ),
where N (γ

(t)
s\K ; Bν (t) , Σ) denotes the multivariate Gaussian distribution, with mean vector Bν (t) and covariance matrix Σ, evaluated at γ

(t) s\K . Finally p(ν|µ 0 , A, Φ, V 0 ) = p(ν (1) |µ 0 , V 0 ) T t=2 log p(ν (t) |ν (t-1) , A, Φ).

Estimation

This section focuses on the inference of the model proposed above. A variational EM algorithm is considered and a model selection criterion is derived.

A variational framework

We aim at maximizing the log-likelihood log p(X|θ) associated with the model. To achieve this maximization, a common approach consists in using an expectation maximization (EM) algorithm (Dempster et al., 1977;Krishnan and McLachlan, 1997). However, such an algorithm cannot be derived here since p(Z, γ, ν|X, θ) is intractable. Therefore, we propose to use a variational EM-type algorithm (VEM) (Hathaway, 1986) which locally optimizes the model parameters with respect to a lower bound of the log-likelihood. Thus, given a distribution q for the three sets of latent variables (Z, γ, ν), the log-likelihood can be written:

log p(X|θ) = L(q, θ) + KL(q(.) p(.|X, θ)), ( 4
)
where L is defined as follows:

L(q, θ) = Z γ ν q(Z, γ, ν) log p(X, Z, γ, ν|θ) q(Z, γ, ν) dγ dν, (5)
and KL denotes the Kullback-Leibler divergence between the true and approximate posterior distributions:

KL(q(.) p(.|X, θ)) = - Z γ ν q(Z, γ, ν) log p(Z, γ, ν|X, θ) q(Z, γ, ν) dγ dν. (6)
Looking for the best approximation of the posterior distribution p(Z, γ, ν|X, θ) in the sense of the KL divergence becomes equivalent to searching for a distribution q(•) that maximizes the lower bound L of the integrated log-likelihood.

Unfortunately, because the joint distribution (2) in the lower bound involves the quantity p(Z|γ) which depends on the normalizing constant C(γ

(t)
s ), L has no analytical form and cannot be optimized with respect to q(•). Indeed, C(γ

(t) s ) = log( K l=1 exp(γ (t) sl )
) is based on a non linear transformation of the vector γ s ∈ R * + , a lower bound of the first lower bound L(q, θ) is given by: log p(X|θ) L(q, θ) L(q, θ, ξ), where L(q, θ, ξ)

= Z γ ν q(Z, γ, ν) log p(X|Z, Π)h(Z, γ, ξ)p(γ \K |B, ν, Σ)p(ν|µ 0 , A, Φ, V 0 ) q(Z, γ, ν) dγdν (7) with log h(Z, γ, ξ) = T t=1 K k=1 N i=1 S s=1 y is Z (t) ik γ (t) sk -ξ -1(t) s K l=1 exp(γ (t) sl )-1+log(ξ (t) s ) .
Note that the variational parameters ξ (t) s can be optimized to obtain tight bounds (see the end of Section 3.2). Moreover, we emphasize that a variational parameter ξ (t) s is considered for each subgraph s and each time t for more flexibility and to improve the inference procedure. We point out that the quality of the variational approximation we propose cannot be tested analytically since L(q, θ, ξ) and the Kullback-Leibler divergence in (6) are not tractable. Nevertheless, we rely on them for inference purposes. Note that similar approximation schemes have been used for instance by Bishop and Svensén (2003) and Latouche et al. (2014), in the context of model selection.

In order to maximize L(q, θ, ξ), we further assume that q(Z, γ, ν) can be factorized:

q(Z, γ, ν) = q(Z)q(γ)q(ν) = T t=1 N i=1 q(Z (t) i ) q(γ)q(ν).
Finally q(γ) is chosen within the family of Gaussian distributions of the form:

q(γ) = T t=1 S s=1 K k=1 N (γ (t) sk ; γ(t) sk , σ2 (t) sk ),
to derive analytical expectations in the E step, as in Lafferty and Blei (2006).

Since the last component of each vector γ (t) s has to remain equal to zero, to preserve the bijectivity constraints of the transformation f (•), the terms γ(t) sK and σ2 (t) sK are all set to zero to ensure a Dirac mass at zero. All other mean and variance terms (γ

(t) sk , σ2 (t)
sk ), ∀s, k = K, t, are parameters to be estimated.

A VEM algorithm for the dRSM model

In this section, we first assume that the variational terms ξ, which were introduced for approximation purposes, are given. This allows the use of a VEM algorithm (Jordan et al., 1999) to maximize the lower bound L(q, θ, ξ) with respect to q(Z, γ, ν) and the model parameters θ. Such an optimization procedure is iterative and involves a series of successive updates. In the E step, the model parameters are fixed and the lower bound is optimized with respect to q(Z, γ, ν). Conversely, during the M step, the variational distribution is held fixed while L(q, θ, ξ) is maximized with respect to θ. In standard VEM algorithms, a unique set of latent variables is usually considered. In our case, there are three sets (Z, γ, ν) of latent variables and therefore the E step itself involves iterative updates (as in Latouche et al., 2014, for instance). All distributions in q(Z, γ, ν) are held fixed, except one, which is optimized. This procedure is repeated for all distributions in turn.

In the following, we give the update formulae for the E and M steps. The details of the calculations along with the derivation of the lower bound are given in the appendix.

Proposition 3.2. The VEM update step for each distribution q(Z (t) i ) is given by: q(Z

(t) i ) ∼ M Z (t) i ; 1, τ (t) i = (τ (t) i1 , . . . , τ (t)
iK ) ∀i, t, where

τ (t) ik ∝ exp K l=1 C c=0 N i =j δ(X (t) ij = c)τ (t) jl log(Π c kl ) + log(Π c lk ) + S s=1 y is γ(t) sk -ξ -1(t) s K l=1 exp(γ (t) sl + σ2 (t) sl 2 ) -1 + log(ξ (t) s ) . Note that τ (t)
ik is the approximate posterior probability that node i belongs to cluster k at time t.

Proposition 3.3. The VEM update step for the distribution q(ν) is given by:

q(ν) ∝ p(ν (1) |µ 0 , V 0 ) T t=2 p(ν (t) |ν (t-1) , A, Φ) T t=1 N ( S s=1 γ(t) s S ; Bν (t) , Σ S ) .
At this step, we recall that the terms γ(t) s are fixed and so is the variable

x (t) = S s=1
γ(t) s /S. Therefore, it is remarkable to note that the functional form of q(ν) corresponds exactly to the form of the posterior distribution associated with a state space model where ν is the set of all latent state variables and x = (x (t) ) t the set of observed outputs. Thus, each x (t) can be written as x (t) = Bν (t) + ṽ where ṽ ∼ N (0, Σ/S) while the variables in ν are defined as previously:

ν (t) = Aν (t-1) + ω ν (1) = µ 0 + u, with ω ∼ N (0, Φ) u ∼ N (0, V 0 ).
Contrary to the original state space model introduced in Section 2, where both γ and ν were sets of unobserved variables, we obtain here a standard linear dynamic system from which the corresponding parameters, i.e. θ = (µ 0 , A, B, Φ, V 0 , Σ/S) can be estimated using Kalman filter and Rauch-Tung-Striebel (RTS) smoother (Rauch et al., 1965) (details can also be found in Minka, 1998). The expectations ν(t) and covariance matrices V (t) of the random variables ν (t) , given all the observed data x, are determined relying on backward forward recursions.

Proposition 3.4. After the E step of the VEM algorithm, the lower bound L(q, θ, ξ) simplifies into:

L(q, θ, ξ) = T t=1 K k,l C c=0 N i =j δ(X (t) ij = c)τ (t) ik τ (t) jl log(Π c kl ) + T t=1 S s=1 r (t) s γ(t) sk -N s ξ -1(t) s K l=1 exp(γ (t) sl + σ2 (t) sl 2 ) + N s -N s log(ξ (t) s ) + T t=1 S s=1 log N (γ (t) s , B ν(t) s , Σ) - 1 2 tr(Σ -1 B V (t) B) - 1 2 tr(Σ -1 σ(t) 2 s ) - T t=1 S s=1 K-1 k=1 -log (2π) 1 2 σ(t) sk + T KS 2 - T t=1 log N (x (t) ; B ν(t) , Σ S ) + 1 2 tr(Σ -1 SB V (t) B) - N i=1 T t=1 K k=1 τ (t)
ik log(τ

(t) ik ) + log p(x|θ )
where r

(t) s = N i=1 τ (t)
ik y is , N s is a number of nodes in the subgraph s, and log p(x|θ ) is the log likelihood of the linear dynamic system associated with the variational distribution q(ν) (see Proposition 3.3)

The maximization of this bound allows to obtain the updating formula for the tensor matrix Π:

Πc kl = T t=1 N i =j δ(X (t) ij = c)τ (t) ik τ (t) jl T t=1 C c=0 N i =j δ(X (t) ij = c)τ (t) ik τ (t) jl , ∀k, l, c.
For the parameters γ(t) sk and σ2 (t) sk , we do not obtain analytical expressions, and therefore we rely on a quasi-Newton algorithm for the optimization task.

Optimization of ξ

So far, we have seen that a VEM algorithm could be implemented from approximations depending on the variational parameters ξ (t) s . However, we have not addressed yet how these parameters could be estimated from the data. We follow the work of Svensén and Bishop (2004) on Bayesian hierarchical mixture of experts. Thus, the lower bound L(q, θ, ξ) is optimized with respect to the variational terms ξ (t) s to obtain the tightest bound L(q, θ, ξ) of L(q, θ). This leads to new estimates ξ(t) s of ξ

(t) s : ξ(t) s = K l=1 exp(γ (t) sl + σ2 (t) sl ), ∀s, t.
This procedure gives rise to a three step optimization scheme. Given all ξ = (ξ

(t)
s ) st , the VEM algorithm described previously is used to maximize the lower bound with respect to q(Z, γ, ν) and θ. These terms are then held fixed and a new estimate of ξ is computed. The three steps are repeated until convergence of the lower bound.

Model selection: choice of the number K of latent groups

Using the VEM algorithm proposed in the previous paragraphs, the estimation of the model parameters and of the group memberships is fully automatic for a given value of K. Since we consider here a model-based approach, two dRSM models with different values of K can be considered as two different models. The problem of choosing K can therefore be viewed as a model selection problem. It can be tackled in a model-based context using model selection criteria, such as the Akaike information criterion (AIC) (Akaike, 1974) or the Bayesian information criterion (BIC) (Schwarz, 1978). Due to its popularity and its asymptotic properties (Leroux, 1992), we use BIC in the numerical experiments presented in the following sections. BIC relies on an asymptotic approximation of the marginal log-likelihood, also called integrated log-likelihood, and is defined in the specific context of the dRSM model M by:

BIC(M) = log p(X| θ) - η(M) 2 log(T N (N -1)),
where η(M) is the number of free model parameters depending on K, for the identifiability constraints considered. Unfortunately, the log-likelihood log p(X| θ) = log Z p(X, Z| θ is not tractable here because it involves marginalizing over all latent vectors Z (t) i in Z. Therefore, we propose to replace the log-likelihood with its variational approximation L(q, θ, ξ). Thus, the VEM algorithm is run for various values of K. For each K, the algorithm iterates until convergence of the lower bound. K is then chosen such that the (approximate) BIC criterion is maximized.

Numerical experiments and comparisons

This section aims at proving on synthetic data the validity of the inference algorithm presented in section 3. An introductory example is first considered to highlight the main features of the proposed approach. Model selection is then considered to validate the criterion choice. Extensive comparisons with state-of-the-art methods conclude this section.

Experimental setup

In order to validate our approach, we use in this section artificial data generated according to a common experimental setup. To simplify the characterization and facilitate the reproducibility of the experiments, we designed five different scenarios. The generation setup for each scenario is summarized in

Parameters Scenario 0 Scenario 1 Scenario 2 Scenario 3 Scenario 4 N 300 K 4 T 10 (indep.) 10 (SSM) S 1 1 1 2 2 C 1 1 1 1 2 (Π 0 ll ) l=1,..,K (0.1,0.4,0.5,0.6) Π 0 kl,k =l 0.99 0.8 0.99 Π c =0 kl (1 -Π 0 kl )/C
Table 2: Parameter values for the five types of graphs used in the experiments. In scenario 0, the networks are drawn without an explicit temporal dependance whereas, in the other scenarios, the temporal dependance is generated through a state space model (SSM).

Table 2. Data from scenario 0 are drawn using SBM at each time t and without an explicit temporal dependence. The data sets for all other scenarios (scenarios 1 to 4) are drawn according to the dRSM model. Therefore, the temporal dependence is generated through a state space model. All generated networks are made of N = 300 nodes, distributed into K = 4 latent groups and have T = 10 time points. Depending on the scenario, the networks have S = 1 or 2 subgraphs, with binary (C = 1) or categorical (C = 2) edges. When S > 1, the nodes are randomly assigned uniformly to the subgraphs. Notice that scenario 2 has a parameter Π 0 kl,k =l equal to 0.8 which leads to less heterogeneous latent groups.

The model parameters used for the simulation are as follows. For the simulation of γ, it is assumed that the matrices A, B and V 0 are set to I K-1 , and that Σ = 0.1 × K × I K-1 and Φ = 0.01 × I K-1 . Finally, the tensor matrix Π, which defines the connection probabilities between clusters for the C different types, is set up such that, within the clusters, the probability 1 -Π 0 ll of having an edge of any type is larger than the corresponding connection probabilities between clusters 1 -Π 0 kl,k =l (see Table 2). Notice that such a choice of parameters induces networks made of communities. Then, in case of a connection between two nodes, the edge type is sampled uniformly, i.e. Π c =0 kl = (1 -Π 0 kl )/C, ∀k, l.

An introductory example

We first focus on an introductory example to illustrate the global behavior of the proposed methodology. To this end, we simulated a single network according to scenario 2 for facilitating the understanding of the results. We remind that in this setup the number K of latent groups is fixed to 4 and that C = 1. Therefore, the network is binary and Π 1 kl indicates the occurrence probability of an edge. We ran the VEM algorithm on it for a number K of groups ranging from 3 to 6. We selected afterward the most appropriate number of groups using the BIC criterion.

Figure 4 shows the BIC values associated to the results provided by our VEM algorithm for the different values of K. One can observe that the criterion picks at K = 4, which is the actual simulated value for K. To quickly assess the estimation quality, Table 3 allows to compare the actual (left panel) and estimated (right panel) values of the terms Π 1 kl in the tensor matrix Π, which define the connection probabilities between the latent clusters. On this single example, the estimated values Π 1 kl turn out to be extremely close to the true ones. Similarly, Figure 6 compares the actual (dashed red lines) and estimated (solid black lines) values of the group proportions α for the simulated example. Once again, the estimation of α appears to be very close to the true proportions.

Choice of K

We now focus on the evaluation of the criterion we proposed to select the number K of latent groups. Since our approach aims at searching the unobserved clustering partition the nodes, we chose here to evaluate the combination of our VEM algorithm with the BIC criterion by comparing the resulting partition with the actual one (the simulated partition). In the clustering community, the adjusted Rand index (ARI) (Rand, 1971) serves as a widely accepted criterion for the difficult task of clustering evaluation. The ARI looks at all pairs of nodes and check wether they are classified in the same group or not in both partitions. As a result, an ARI value close to 1 means that the partitions are similar and, in our case, that the VEM algorithm succeeds in recovering the simulated partition.

To validate the combination of our VEM algorithm with the BIC criterion, the analysis was repeated for 50 different data sets, generated according to scenario 2, for a number K of latent groups ranging from 3 to 6. This allows us q q q q q q q q q q 2 4 6 8 10 -197970 -197950

Evolution of the bound for K=4

Iterations

Value of the bound to both verify the consistency of the BIC criterion and to study the clustering ability of our approach. Figure 7 shows the repartition of the criterion values (left panel) as well as the associated ARI values (right panel). These results first confirm that BIC is a valid criterion for selecting the number of groups in this context. Indeed, the value K = 4 is the one which is the most frequently associated with the highest value of BIC. We remind that K = 4 is the actual number of latent groups. One can also observe that the partition resulting from our VEM algorithm is associated, for this value of K, to an ARI value extremely close to 1 which denotes a good matching with the actual partition of the data.

Comparison with the other stochastic models

Our third set of experiments now aims at comparing the performance of our approach to that of state-of-the-art methods. We are here interested in the comparison of dRSM with the following methods: SBM (Nowicki and Snijders, 2001), RSM (Jernite et al., 2014) anddM3SBM (Ho et al., 2011). Once again, the evaluation of the results is done using the ARI criterion. In order to fit a SBM on a dynamic network, we ran the mixer package (Ambroise et al., 2010) for the R software at each time t and the ARI is then computed on the concatenation of all group labels. However, let us notice that SBM was not able to handle networks with categorical edges (scenario 3). For RSM, we used the Rambo package (Bouveyron et al., 2013) for R, on an aggregated version of the whole network. Conversely to SBM, RSM is only able to deal with categorical networks and, consequently, it works only in scenario 4. Finally, we used the Matlab toolbox dM3SBM, kindly provided by the authors, to fit the dM3SBM on the dynamic networks. However, dM3SBM is also not able to handle networks with categorical edges (scenario 4). In order to consider a wide type of networks, we compare here the methods over the five simulation scenarios. We remind that Table 2 summarizes the main features of each scenario. This comparison has been conducted in two different situations: with and without the knowledge of the actual number of clusters. Table 4 presents the clustering results for the four studied methods in the case where the actual number K = 4 of groups has been provided to each method. Conversely, Table 5 presents the clustering results when the methods have to look for the value of K. Reported values are averaged ARI values (with standard deviations) on 20 networks for each scenario. The average selected number K of latent groups is also provided for Table 5.

First, for scenarios 0, 1 and 2, which consider dynamic networks with binary edges (C = 1) and with only one subgraph (S = 1), one can see on Tables 4 and5 that SBM is, as expected, not able to handle the network dynamic. Indeed, SBM obtains a low ARI value in all situations, even though it correctly estimates the number of clusters (Table 5). Conversely, the two dynamic methods (dM3SBM and dRSM) turn out to be able to recover the clustering structure of the dynamic networks. One can however notice that dRSM significantly outperforms dM3SBM in this situation. Notice also the accurate estimation of the number K of clusters made by dRSM (Table 5).

In scenario 3, the simulated dynamic networks are now made of two subgraphs (S = 2), still with binary edges (C = 1). Naturally, SBM does not perform well in this situation too. The dM3SBM provides clustering results similar to the ones of previous scenarios: it globally succeeds in recovering the dynamic but fails in recognizing the clustering pattern. On the other hand, dRSM provides again accurate clustering results associated with good estimations of K, meaning that it succeeds in identifying both the dynamic and clustering patterns.

Finally, scenario 4 considers the case of dynamic networks with two subgraphs (S = 2) and categorical edges (C = 2). Only RSM and dRSM are able to deal with this kind of networks. Similarly to SBM in previous scenarios, RSM does not succeed in recovering the dynamic and provides very unsatisfactory clustering results. Conversely, dRSM gives very good clustering results regarding the difficulty of the situation. It is worth noticing the sharp estimation made by dRSM of the number K of group in this case too. This confirms the efficiency of both our inference algorithm and our model selection criterion.

We also used scenario 4 to highlight that providing the methodology with the right subgraph structure helps in clustering the vertices. Thus, with the knowledge of the actual number of clusters, we ran dRSM with the wrong subgraph structure (S = 1), and we obtained an average ARI of 0.54±0.2. This result is to be compared to the ARI performances for scenario 4, as presented in Table 4.

Maritime network

This section presents an application of the proposed methodology for the analysis of a network of maritime flows in which a temporal dynamic is present. The dynamic network was provided by Dr. César Ducruet, from the Géographie-Cités laboratory, who is interested in studying the evolution of maritime flows over time. The data was extracted from the well-known Lloyd's list which has recorded almost all ship movements worldwide since 1890.

Data and study protocol

Data was obtained from the printed Lloyds voyage record published every October between 1890 until 2008. The list details, for each merchant vessel, its successive movements from one port to another. From the raw database of vessel flows, we extracted a dynamic network with 17 time points. The first observation is October 1890 and the network ends in October 2008. Table 6 provides the correspondence between the 17 time points and the actual dates.

At each time point, the adjacency matrix between ports was constructed as follows. First, for every pair of ports, we calculated the total number of ship movements between those ports. Then, we set the associated entry in the adjacency matrix to 1 if the number of ship movements between the two ports is greater or equal to 1, and to 0 otherwise. The original network contained 4472 ports worldwide. We however had to reduce the network size to only 286 ports since most of the ports were not active throughout the whole period of the study. We finally applied dRSM to a maritime network which describes the navigation of ships among 286 ports in the world at 17 time points. Let us highlight that the study period includes many major historical or economical events (the two world wars, the oil crisis, the economic crisis in Europe, ...), which could directly affect the navigation movements at a global scale and could also change the port behaviors.

The partition of the network into subgraphs is here provided by the port memberships to the four main maritime basins: Asia -Pacific, Europe -Atlantic, Mediterranean -Black Seas, and Middle East -Indian Ocean. Figure 8 presents this partition of the ports where the colors indicates the different subgraphs.

To summarize, the network is a undirected and binary network without self loops, i.e. C = 1 and X t ij = 1 if the port i and the port j exchange at least one ship during the period t, 0 otherwise, with t ∈ {1, . . . , 17} and S = 4. Figure 9 shows the adjacency matrix, in 1890 and in 2008, between the 286 ports organized by subgraph.

Results

We used the variational EM algorithm introduced in the previous section in order to find the latent groups that may be hidden in the data. The choice of the number of groups is made by applying the VEM algorithm for K = 3, ..., 8 and by then computing the associated BIC values. The retained value for K is the one associated with the highest BIC value. To ensure a good accuracy of the results, the VEM algorithm was run 5 times for each value of K. Figure 10 shows the evolution of the BIC criterion according to K. One can observe that BIC peaks at K = 7, meaning that 7 latent groups seem to organize the network. We therefore chose this specific value for K and retained the best run for K = 7 over the five runs as the final clustering result.

Network in 1890

Network in 2008 On the one hand, it is of main interest to look at the estimated tensor matrix Π in order to understand and characterize the found latent groups. Indeed, the tensor matrix Π describes the connection probabilities between the groups and allows to figure out the different connection patterns. Since the network considered here is binary, it is enough to look at the terms Π 1 kl since Π 0 kl + Π 1 kl = 1, for all k, l. Figure 11 presents those estimated values. From the figure, clusters 6 and 7 appear to be groups of hubs for which the connection probabilities are large within and between clusters.

On the other hand, the estimated group proportions over time should allow to understand the dynamic of the network. Figure 12 presents the evolution of those proportions over time for each subgraph. One can first observe that the proportion of cluster 6 is low and rather stable over time. This confirms that cluster 6 is a group of a limited number of hubs with a high connectivity and probably a high level of traffic. Cluster 6 includes ports such as Anvers, Rotterdam or Singapore. It is also interesting to see that, in subgraph 2 (Europe -Atlantic), the number of hubs increased until 1930, was then perturbed during the second world war and finally decreased from 1951. Conversely, in subgraph 1 (Asia -Pacific), the proportion of hubs was low until 1975 and then significantly increased. From a global point of view, one can also observe a clear and recent reorganization of the network in which hubs tend to be less numerous worldwide (and probably bigger).

Regarding cluster 7, one can see on Figure 12 that its proportions in the subgraphs are higher than those of cluster 6. The ports of cluster 7 can be qualified as hubs of second class which are subordinated to the main hubs of cluster 6. Most of them are marked by a colonial logic, such as Marseille, Connection probabilities between clusters Kolkata or Cape Town. The evolution of this cluster until the recent period shows a persisting link North-South (e.g. Le Havre -Casablanca) or East-West (e.g. Spain -Brazil -Canaries).

The cluster 5 is mainly made of ports from the Asia -Pacific and Middle East -India basins except during major crises, such as World War II and the oil crisis. During those crises, the cluster mainly contains European ports. The rapid modification of this cluster appears clearly on Figure 12 around 1946Figure 12 around , 1980Figure 12 around and 2008. This cluster can be interpreted as made of active ports from the developing world which move to cluster 2 during the crises. This may highlight the disintegration of long distance links during such crises. Conversely, cluster 2 turns out to be mostly made, except during crises, of European ports of average size, mainly on the atlantic coast. Those ports are rather a reflection of a past glory and most of them have declined over the century. This may due to a failed industrialization or a significant distance to the major trade routes.

Finally, clusters 3 and 4 are made of very small ports with low activity. Those ports are usually not connected together and communicate with the rest of the network only through ports of clusters 2 and 5. The connection with clusters 2 and 5 explains the brutal changes in the proportions of clusters 3 and 4 that one can also observe.

Conclusion

This work has considered the problem of analyzing dynamic networks with categorical edges and for which a subgraph partition is known. This kind of networks is frequent in a wide range of scientific fields, such as Geography in particular. For this purpose, we proposed an extension of the RSM model to the dynamic setting. The new model, called dRSM, uses a state space model to model the evolution of the latent group proportions over time. A variational expectation maximization (VEM) algorithm is proposed to perform inference. We have shown in particular that the variational approximations lead to a new state space model from which the parameters can be estimated using the standard Kalman filter and the Rauch-Tung-Striebel (RTS) smoother. Model selection is also considered through an approximate BIC criterion.

Numerical experiments have highlighted the main features of the dRSM model and have demonstrated the efficiency of both the VEM algorithm and the model selection criterion. A numerical comparison has also shown that existing methods, dynamic or not, are less flexible and efficient than dRSM when applied to dynamic networks. Finally, dRSM has been applied to a dynamic maritime flow network, build from the famous Lloyd's list, and has allowed to characterize interesting dynamic phenomena. painstaking analysis of the results. The data were collected in the context of the ERC Grant N.313847 "World Seastems" (http://www.world-seastems.cnrs.fr). The authors would like also to thank Catherine Matias and Stéphane Robin for their useful remarks and comments on this work.

A. Construction of a tractable lower bound

We rely on a bound introduced in Jordan et al. (1999). Such a general bound can easily been derived by noticing that C(•) is a concave function of

K l=1 exp(γ (t)
sil ) and therefore a first order Taylor expansion of the normalizing constant, at any ξ (t) s ∈ R * + , will lead to the inequality:

log( K l=1 exp(γ (t) sl )) ≤ ξ -1(t) s ( K l=1 exp(γ (t) sl )) -1 + log(ξ (t) s ). (8)
The bounds (8) on the C(γ

(t)
s ) terms induce a lower bound on the quantity log p(Z|γ):

log p(Z|γ) = T t=1 K k=1 N i=1 Z (t) ik log(f k (γ (t) si )) = T t=1 K k=1 N i=1 S s=1 y is Z (t) ik γ (t) sk -log( K l=1 exp(γ (t) sl )) ≥ log h(Z, γ, ξ),
where ξ denotes the set of all variational parameters (ξ

(t) s ) st and the function h(•, •, •) is such that: log h(Z, γ, ξ) = T t=1 K k=1 N i=1 S s=1 y is Z (t) ik γ (t) sk -ξ -1(t) s K l=1 exp(γ (t) sl )-1+log(ξ (t) s ) .
Replacing log p(Z|γ) by log h(Z, γ, ξ) in L(q, θ), leads to a new lower bound L(q, θ, ξ) for log p(X|θ) which satisfies: log p(X|θ) L(q, θ) L(q, θ, ξ), where L(q, θ, ξ)

= Z γ ν q(Z, γ, ν) log p(X|Z, Π)h(Z, γ, ξ)p(γ|B, ν, Σ)p(ν|µ 0 , A, Φ, V 0 ) q(Z, γ, ν) dγdν.

B. E-step of the VEM algorithm

Distribution q(Z)

The VEM update step for each of the distributions q(Z i ) in q(Z) is given by: log

q(Z i ) = E γ,ν,Z \i [log p(X|Z, Π) + log h(Z, γ, ξ)] + const = T t=1 K k=1 Z (t) ik   K l=1 C c=0 N i =j δ(X (t) ij = c)τ (t) jl log(Π c kl ) + log(Π c lk )   + T t=1 K k=1 S s=1 y is E γ Z (t) ik log h(Z (t) , γ (t) , ξ (t) ) + const. = T t=1 K k=1 Z (t) ik   K l=1 C c=0 N i =j δ(X (t) ij = c)τ (t) jl log(Π c kl ) + log(Π c lk )   + T t=1 K k=1 S s=1 y is E γ γ (t) sk -(ξ -1(t) s K l=1 exp(γ (t) sk ) -1 + log(ξ (t) s )) + const. = T t=1 K k=1 Z (t) ik   K l=1 C c=0 N i =j δ(X (t) ij = c)τ (t) jl log(Π c kl ) + log(Π c lk )   + T t=1 K k=1 S s=1 Z (t) ik y is γ(t) sk -ξ -1(t) s K l=1 E(exp(γ (t) sl )) -1 + log(ξ (t) s ) + const. = T t=1 K k=1 Z (t) ik K l=1 C c=0 N i =j δ(X (t) ij = c)τ (t) jl log(Π c kl ) + log(Π c lk ) + S s=1 y is γ(t) sk -ξ -1(t) s K l=1 exp(γ (t) sl + σ2 (t) sl 2 ) -1 + log(ξ (t) s ) + const,
where all terms that do not depend on Z i have been put into the constant terms const. Moreover since γ

(t) sk ∼ N (γ (t)
sk , σ2 (t) sk ) we have used:

E[exp(γ (t)
sk )] = exp(γ

(t) sk + σ2 (t) sk 2 ).
We then recognize the functional form of a multinomial distribution:

q(Z (t) i ) ∼ M(Z (t) i ; 1, τ (t)
i ), ∀i, t.

Distribution q(ν):

The VEM update step for the distribution q(ν) is given by: log q(ν) = E Z,γ log p(γ|ν, Σ, B) + log p(ν|µ 0 , V 0 , A, Φ) + const

E γ - 1 2 (γ (t) s ) Σ -1 (γ (t) s ) + (γ (t) s ) Σ -1 Bν (t) - 1 2 (ν (t) ) B Σ -1 Bν (t)
+ log p(ν (1) |µ 0 , V 0 ) + + log p(ν (1) |µ 0 , V 0 ) + T t=2 log p(ν (t) |ν (t-1) , A, Φ) + const, where all terms that do not depend on ν have been put into the constant terms const. We recognize the functional form of the posterior distribution of a linear dynamic system: log q(ν) = log p(ν (t) |ν (t-1) , A, Φ) + const.

C. Derivation of the lower bound

In the following, we denote x (t) = S s=1 γ(t) s S an observed variable. The lower bound is given by: L(q, θ, ξ) = Z γ ν q(Z, γ, ν) log p(X|Z, Π)h(Z, γ, ξ)p(γ|ν, Σ)p(ν|µ 0 , A, Φ, V 0 ) q(Z, γ, ν) dνdγ = E Z,γ,ν [log p(X|Z, Π)h(Z, γ, ξ)p(γ|ν, Σ, B)p(ν|µ 0 , A, Φ, V 0 ) q(γ)q(ν) N i=1 q(Z i ) ] = E Z (log p(X|Z, Π)) + E Z,γ (log h(Z, γ, ξ)) + E γ,ν (log p(γ|ν, Σ, B))

+ E ν (log p(ν|µ 0 , A, Φ, V 0 )) -E γ (log q(γ)) -E ν (log q(ν)) -E Z (log( N i=1 q(Z i ))).

Note that (see Proposition 3.3), q(ν) ∝ p(ν (1) |µ 0 , V 0 ) T t=2 p(ν (t) |ν (t-1) , A, Φ) As pointed out in this proposition, this corresponds to the form of the posterior distribution associated with a state space model with parameter θ and with observed outputs x = (x (t) ) t . If we denote p(x|θ ) the likelihood associated with this model, and the joint likelihood p(x, ν|θ ), we have q(ν) = p(x, ν|θ ) p(x|θ ) .

Therefore E ν (log q(ν)) = E ν (log p(ν|µ 0 , A, Φ, V 0 ))+E ν (log p( and L(q, θ, ξ) can be written as follows:

L(q, θ, ξ) = Z γ ν q(Z, γ, ν) log p(X|Z, Π)h(Z, γ, ξ)p(γ|ν, Σ)p(ν|µ 0 , A, Φ, V 0 ) q(Z, γ, ν)

= E Z,γ,ν [log p(X|Z, Π)h(Z, γ, ξ)p(γ|ν, Σ, B)p(ν|µ 0 , A, Φ, V 0 ) q(γ)q(ν) q(Z i )))

+ log p(x|θ ).

We explicit below each of the terms of the bound L(q, θ). log N (x (t) ; B ν(t) , Σ/S) -1 2 tr(Σ -1 SB T V (t) B) .

6. E Z (log( T i=1 q(Z i ))):

Figure 1 :

 1 Figure 1: Connections between a subset of 26 ports (from October 1890 to October 2008). Data extracted from Lloyd's list. The known subgraphs correspond to geographical regions (continents) indicated using colors.

Figure 3 :

 3 Figure 3: Graphical representation of the dRSM model.

  (t) swhich makes some expectations of the standard VEM algorithm impossible to derive.Following the work of Lafferty and Blei (2006) on correlated topic models, we propose a new bound of L(q, θ) based on a variational lower bound of p(Z|γ), as inJordan et al. in Jordan et al. (1999). Proposition 3.1. (Proof in A) Given any set ξ of variational parameters ξ (t)

FigureFigure 4 :

 4 Figure 4: Choice of K by model selection with BIC for a simulated network. The actual value for K is 4.

Figure 5 :

 5 Figure 5: Evolution of the bound L for K = 4.

Figure 7 :

 7 Figure 7: Criterion and ARI values over 50 networks generated.

Figure 8 :

 8 Figure 8: The given partition of the 286 nodes (ports) into 4 subgraphs.

Figure 9 :Figure 10 :

 910 Figure 9: Adjacency matrix of the maritime network organized by subgraph (basin) in 1890 (left) and 2008 (right).

Figure 11 :

 11 Figure 11: Terms Π 1 kl of the tensor matrix Π estimated using the VEM algorithm.

E

  γ log N (γ (t) s ; Bν (t) , Σ) + log p(ν (1) |µ 0 , V 0 )

  ν (t) |ν (t-1) , A, Φ) + const. t) ) B (SΣ -1 )Bν (t)

  p(ν (1) |µ 0 , V 0 ) + T t=2

  )-log p(x|θ ).This leads to,E ν (log p(ν|µ 0 , A, Φ, V 0 ))-E ν (log q(ν)) = -E ν (log p(

  (log p(X|Z, Π)) + E Z,γ (log h(Z, γ, ξ)) + E γ,ν (log p(γ|ν, Σ, B)) -E γ (log q(γ)) -E ν (log p(

  1. E Z (log p(X|Z, Π)):E Z (log p(X|Z, Π)) = Z,γ (log h(Z, γ, ξ)): E Z,γ (log h(Z, γ, ξ)) = E Z

  )

	Notations Description
	X	Adjacency matrix X ij ∈ {0, . . . , C} at each t (t)
	Z	Binary matrix. Z	(t) ik = 1 indicates that i belongs to cluster k at t
	N	Number of vertices in the network
	K	Number of latent clusters
	S	Number of subgraphs
	C	Number of edge types
	Π	Π c kl is the probability of having an edge of type c between
		vertices of clusters k and l
	α	α

Table 4 :

 4 Clustering results for the four studied methods on networks simulated according to the five scenarios. The actual number K = 4 of groups has been provided to each method here. Average ARI values are reported (with standard deviations) and results are averaged on 20 networks for each scenario.

	Method	Scenario 0	Scenario 1	Scenario 2	Scenario 3 Scenario 4
	SBM	0.10±0.04	0.12±0.05	0.18±0.07	0.14±0.09	-
	RSM	-	-	-	-	0.01±0.01
	dM3SBM	0.36±0.09	0.30±0.16	0.25±0.16	0.32±0.20	-
	dRSM	1.00±0.00	0.98±0.04	0.90±0.20	0.97±0.07	0.75±0.24

Table 6 :

 6 The time points considered in the maritime network.

  E γ,ν (log p(γ|ν, Σ, B)): E γ,ν (log p(γ|ν, Σ, B)) = E γ,ν log

	3. T			S
											N (γ (t) s ; Bν (t) s , Σ)
							t=1	s=1
			=	T t=1	S s=1		log N (γ (t) s , B ν(t) s , Σ) -	1 2	tr(Σ -1 B T V (t) B) -	1 2	tr(Σ -1 σ(t) 2 s )
	4. E γ (log q(γ)):
					T		S			K
	E γ (log q(γ)) = E γ							N (γ	(t) sk ;	γ(t) sk , σ2 (t) sk )
					t=1	s=1	k=1
			=	T t=1	S s=1	K k=1	-log (2π)	1 2 σ(t) sk -	T KS 2	.
	5. E ν (log p(	S s=1	γ(t)
		S s=1	γ(t)
											y is τ ik (t)	γ(t) sk -τ ik ξ -1(t) (t) s	K l=1	exp(γ	(t) sl +	σ2 (t) sl 2	)
			+ τ ik -τ (t) ik log(ξ (t) (t) s )
			=	T t=1	S s=1		r (t) s	γ(t) sk -N s ξ -1(t) s	K l=1	exp(γ	(t) sl +	σ2 (t) sl 2	) + N s -N s log(ξ (t) s )
	where denote r s is a quantity (t)	N i=1 τ	(t) ik y is .

s S |ν (t) , Σ S , B)) : E ν (log p( s S |ν (t) , Σ S , B)) = T t=1

Acknowledgments

The authors would like to greatly thank César Ducruet, from the Géographie-Cités laboratory, Paris, France, for providing the maritime network and for his


Florian Jeanne
email: florian.jeanne@hds.utc.fr

Yann Soullard
email: yann.soullard@hds.utc.fr

Indira Thouvenin
email: indira.thouvenin@hds.utc.fr

What is wrong with your gesture ? An error-based assistance for gesture training in virtual environments

Keywords: 3D interaction, Immersion, Training, Gesture

Virtual environments are commonly used to train people, especially for gesture training. A major challenge consists in error perception. In this paper, we propose to guide the user gesture with an errorbased feedback in a 3D virtual environment. In a preliminary experiment, trainees reproduce an expert gesture from an aluminium foundry. Their gesture performance is evaluated according to two critical error criteria to improve the interaction for training.

INTRODUCTION

In an industrial context, new technicians are generally trained by an expert, who first presents the tasks and gestures of the technical process and then assists them in the first stages of training. Such training is performed on a real workstation in the factory and may lead to safety problems due to trainees lack of expertise. Besides, they may cause some damage on apparatus, leading to important costs for the company.

On the last few years, gesture training systems using virtual environments have been proposed for specific applications such as medical science [START_REF] Moreau | A method to evaluate skill transfer and acquisition of obstetric gestures based on the curvatures analysis of the position and the orientation[END_REF][START_REF] Papadopoulos | Practical chirocentric 3dui platform for immersive environments[END_REF], non-verbal communications [START_REF] Ravenet | From a user-created corpus of virtual agent non-verbal behavior to a computational model of interpersonal attitudes[END_REF] or technical gestures [START_REF] Da Dalto | The use and benefits of virtual reality tools for the welding training[END_REF]. However, most of systems do not provide any kind of feedback regarding errors made by a trainee.

To improve industrial gesture training, we propose a training system based on learning through error with a virtual environment. The gesture error is measured using a machine learning method. Then, sensory feedback are generated in a 3D virtual environment. The machine learning method is a Dynamic Time Warping algorithm (DTW) [START_REF] Berndt | Using dynamic time warping to find patterns in time series[END_REF] which enables to evaluate the matching between a trainee gesture and an expert one on error criteria. The DTW computes a distance for each error criterion (e.g. speed, orientation): the more important the distance is, the worst the error is. Each error criterion refers to a particular 3D visual feedback.

We illustrate our approach with a specific technical gesture from an assembly line in an aluminium foundry. Trainees have to reproduce an expert gesture for mold cleaning with a compressed-air blower. We observe two types of errors: gestures speeds and blower inclination errors. Gestures are tracked using infrared cameras detecting marker positions and orientations.

TRAINING SYSTEM IN VIRTUAL ENVIRONMENT

We identify the most important error to be corrected and then provide an appropriate guidance to improve trainee gestures.

Gesture analysis based on error recognition

The motion tracking system we use extracts marker positions and orientations. However, these features may not be adequate to recognize different types of errors. Thus, we compute higher level features refering to spatio-temporal descriptors and allowing a better recognition of the different error classes.

Working with technical gestures leads to various constraints limiting the choice of a model. Indeed, one generally has to deal with only a few number of gesture examples since they must to be manually labeled by an expert, which is costly. Besides, one wants to get interpretable outputs in order to provide relevant feedback to trainees. Based on these comments, a DTW algorithm [START_REF] Berndt | Using dynamic time warping to find patterns in time series[END_REF] seems to be an appropriate choice for such a task. Indeed, a DTW may be efficient on small training samples, in particular with one gesture template. In addition, it provides interpretable results in the form of a distance between two time series. Besides, unlike many learning methods, a DTW directly works with labeled examples and it does not model the data distribution (which creates a bias).

A DTW is defined to find an optimal alignment between two time series. It computes a distance between the two series based on the optimal alignment. Let x = (x 1 , ..., x N ) and q = (q 1 , ..., q M ) be two time series. The DTW algorithm consists in computing a cumulative cost C(n, m) such that C(n, 1) = ∑ n i=1 d(x i , q 1 ), C(1, m) = ∑ m i=1 d(x 1 , q i ) and for 1 < n N and 1 < m M the cost is computed using the following equation:

C(n, m) = min{C(n-1, m),C(n-1, m-1),C(n, m-1)}+d(x n , q m )
(1) where d(x n , q m ) measures the distance between the n th value of x and the m th value of q. Initial and final conditions are defined to ensure that the process begins and ends with the alignment of the first and last elements of q and x.

For each error criterion, we have a unidimensional time series refering to the feature values computed from the marker positions and orientations at each frame. We apply one DTW per error criterion (speed and inclination): each of them evaluates the matching between a trainee gesture and an expert gesture and provides a distance which informs on the gap between the two gestures for the given type of error. These measures are then provided to the assistance module in order to select an appropriate guidance.

Error-based assistance

The guidance can be multimodal and be related to different nature of feedback such as visual [START_REF] Moreau | A method to evaluate skill transfer and acquisition of obstetric gestures based on the curvatures analysis of the position and the orientation[END_REF], auditory [START_REF] Boyer | Legos project -state of the art[END_REF] or haptic [START_REF] Lindeman | The tactapack: A wireless sensor/actuator package for physical therapy applications[END_REF]. Besides, Galligan et al. [START_REF] Galligan | Advanced PE for Edexcel[END_REF] define different types of feedback: there are intrinsic feedback from the kinaesthetic senses, extrinsic feedback such as performance or results, positive and negative feedback to inform users on the performance level, terminal feedback which are provided before or after the performance, and finally concurrent feedback which are given during the performance. The authors point out that the type of feedback may depend on trainees themselves or on gestures they have to learn, e.g. some gestures require real-time feedback whereas others require an a posteriori feedback.

At the end of a trainee performance, the higher distance provided by the error recognition module refers to the most critical error. Thus, the guidance refering to this most critical error is given to the user. We provide extrinsic feedback enabling users to visualize their performance. With such feedback, we expect trainees to improve their gestures by observing their performance. When the speed is the most critical error, we display the gesture trajectory using spheres whereas we use capsules for an inclination error (Figure 1). Capsules are orientated in the same direction as the tracked marker during the gesture. A color coding informs trainees on the differences between their gesture and the reference one. This is determined according to the optimal alignment provided by the DTW: sky blue and dark blue refer respectively to a slightly and much smaller value than the expert gesture value, green is for a close value whereas orange and red refer respectively to a slightly and much higher value.

PRELIMINARY EXPERIMENT

We applied our approach to the context of aluminium foundry. Within this context, technicians perform various technical gestures. Also, as a first step, we selected one of these gesture which consist of cleaning a mold using a compressed-air blower.

Experimental design

The task consists in cleaning 4 lines of 5 dirts, as illustrated in Figure 1). An expert gesture is beforehand recorded. It consists in cleaning the mold using the blower from the front to the back and from the left to the right with a constant speed when the blower moves from the front to the back of the mold. It has to be performed at a proper speed, otherwise technicians would waste time if they are too slow or the mold would not be correctly cleaned if they are too quick. Finally, the blower has to be lifted up when the technician reaches the back of the mold.

We evaluate users capability to reproduce the technical gesture performed by the expert. Thus, we use a multidimensional DTW which measures the matching between two multidimensional time series. Here gestures are 2-dimensional sequential data where the two dimensions refer to the two error criteria (the speed and the inclination of the blower).

Experimental procedure

The expert begins the process by describing the gesture to the trainee. Then users have to reproduce it 10 times. Between each training iteration, visual feedback (one of the two illustrated in Figure 1) is provided to the trainee according to the error criterion having the higher distance with the expert gesture (i.e. speed or inclination).

Participants and platform

15 users contributed to our experiment, 12 men and 3 women from 18 to 53. All the participants were novices with our training system. We used 6 infrared cameras: 2 OptiTrack Prime 13W and 4 Op-tiTrack Prime 13 (optitrack.com), both with a frequency of 240 Hz ensuring a low latency. A stereoscopic display of a mold with dirts were done on a horizontally laid screen. A compressed-air blower enables to clean the mold remotely. Markers were fixed at the end of the prop and on the head and they were tracked using the system. User performance generally improves during the training process. A Wilcoxon test between the first and the last training iteration shows a significant improvement (the p-value is equal to 0.038 which is significant to a standard 0.05 threshold). In addition, the standard deviation slightly decreases. This illustrates that the gap between user performance decreases during the training. Finally, we point out that user performances seem to remain constant on the last iterations. This may be due to a loss of concentration after some training iterations.

RESULTS AND DISCUSSION

To conclude, those first results seem promising and show that trainee gestures get closer to the expert gesture as expected. In our future works, we will experiment our system on more complex gestures and compared it to other automatic strategy of feedback. Besides, we will focus our research on adaptive dynamic guidances to provide customized assistances.

Figure 1 :

 1 Figure 1: On the left: Illustration of a gesture being performed (red arrows show the directions). On the right: two types of feedback provided a posteriori (the speed criteria represented by balls and the inclination one by capsules)

Figure 2 :

 2 Figure 2: Average gap (distance computed by the DTW) of user performances with the expert one during training iterations.

ACKNOWLEDGEMENTS

T his work is founded by the European Union and the French region Picardie. Europe acts in Picardie with the European Regional Development Fund (ERDF).