Antonino Freno
email: antonino.freno@inria.fr

Mikaela Keller
email: mikaela.keller@inria.fr

Gemma C Garriga
email: gemma.garriga@inria.fr

Marc Tommasi
email: marc.tommasi@inria.fr

Spectral Estimation of Conditional Random Graph Models for Large-Scale Network Data

Generative models for graphs have been typically committed to strong prior assumptions concerning the form of the modeled distributions. Moreover, the vast majority of currently available models are either only suitable for characterizing some particular network properties (such as degree distribution or clustering coefficient), or they are aimed at estimating joint probability distributions, which is often intractable in large-scale networks. In this paper, we first propose a novel network statistic, based on the Laplacian spectrum of graphs, which allows to dispense with any parametric assumption concerning the modeled network properties. Second, we use the defined statistic to develop the Fiedler random graph model, switching the focus from the estimation of joint probability distributions to a more tractable conditional estimation setting. After analyzing the dependence structure characterizing Fiedler random graphs, we evaluate them experimentally in edge prediction over several real-world networks, showing that they allow to reach a much higher prediction accuracy than various alternative statistical models.

INTRODUCTION

Arising from domains as diverse as bioinformatics and web mining, large-scale data exhibiting network structure are becoming increasingly available. Network models are commonly used to represent the relations among data units and their structural interactions. Recent studies, especially targeted at social network modeling, have focused on random graph models of those networks. In the simplest form, a social network is a configuration of binary random variables X uv such that the value of X uv stands for the presence or absence of a link between nodes u and v in the network. The general idea underlying random graph modeling is that network configurations are generated by a stochastic process governed by specific probability laws, so that different models correspond to different families of distributions over graphs.

The simplest random graph model is the Erdős-Rényi (ER) model [START_REF] Erdős | On Random Graphs, I[END_REF], which assumes that the probability of observing a link between two nodes in a given graph is constant for any pair of nodes in that graph, and it is independent of which other edges are being observed. In preferential attachment models [START_REF] Barabási | Emergence of scaling in random networks[END_REF], the probability of linking to any specified node is proportional to the degree of the node in the graph, leading to "rich get richer" effects. Small-world models [START_REF] Watts | Collective dynamics of 'small-world' networks[END_REF]] try to capture instead some phenomena often observed in real networks such as high clustering coefficients and small diameters [START_REF] Mark | Networks. An Introduction[END_REF]. A sophisticated attempt to model complex dependences between edges in the form of Gibbs-Boltzmann distributions is made by exponential random graph (ERG) models [START_REF] Tom | New Specifications for Exponential Random Graph Models[END_REF], which subsume the ER model as a special case. Finally, a recent attempt at modeling real networks through a stochastic generative process is made by Kronecker graphs [START_REF] Leskovec | Kronecker graphs: An approach to modeling networks[END_REF], which try to capture phenomena such as heavy-tailed degree distributions and shrinking diameter properties while paying attention to the temporal dynamics of network growth.

While some of these models behave better than others in terms of computational tractability, one basic limitation affecting all of them is what we may call a parametric assumption concerning the probability laws underlying the observed network properties. In other words, currently available models of network structure assume that the probability distribution gener-ating the network can be expressed through a particular parametric form P (G| θ), where G is the observed graph and θ is a parameter vector. For example, typical formulations of exponential random graph models assume that the building blocks of real networks are given by such structures as k-stars and k-triangles, with different weights assigned to different structures, whereas Kronecker graphs assume that the number of edges and the number of nodes in a given network are related by a densification power law with a suitable densification parameter. In such frameworks, estimating the model from data reduces to fitting the model parameters, whereas the model structure remains fixed from the very beginning. The problem is that, in order for a parametric model to deliver an accurate estimate of the distribution at hand, its prior assumption concerning the form of the modeled distribution must be satisfied by the given data, which is something that we can rarely assess a priori. To date, the knowledge we have concerning real-world network phenomena does not allow to assume that any particular parametric assumption is really capturing in depth any network-generating law, although some observed properties may happen to be modeled fairly well.

The aim of this paper is twofold. On the one hand, we take a first step toward nonparametric modeling of random networks by developing a novel network statistic, which we call the Fiedler delta statistic. The Fiedler delta function allows to model different graph properties at once in an extremely compact form. This statistic is based on the spectral analysis of the graph Laplacian. In particular, it is based on the smallest non-zero eigenvalue of the Laplacian matrix, which is known as Fiedler value [START_REF] Fiedler | Algebraic connectivity of graphs[END_REF][START_REF] Mohar | The Laplacian Spectrum of Graphs[END_REF]. On the other hand, we systematically adopt a conditional approach to random graph modeling, i.e. we focus on the conditional distribution of edges given some neighboring portion of the network, while setting aside the problem of estimating joint distributions. The resulting conditional random graph model is what we call Fiedler random graph (FRG). Roughly speaking, to model the conditional distribution of an edge variable X uv with respect to its neighborhood, we compute the difference in Fiedler values for the vicinity subgraph including or excluding the edge {u, v}. The underlying intuition is that the variations encapsulated in the Fiedler delta for each particular edge will give a measure of the role of that edge in determining the algebraic connectivity of its neighborhood. As part of our contributions, we theoretically prove that FRGs capture edge dependencies at any distance within a given neighborhood, hence defining a fairly general class of conditional probability distributions over graphs.

Experiments on the estimation of (conditional) link distributions in large-scale networks show that FRGs are well suited for estimation problems on very large networks, especially small-world and scale-free networks. Our results reveal that the FRG model supersedes other approaches in terms of edge prediction accuracy, while allowing for efficient computation in the large-scale setting. In particular, it is known that the computation of the Fiedler delta scales polynomially in the size of the analyzed neighborhood [START_REF] Bai | Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide[END_REF]. And the experiments show that even for small sizes of neighborhoods (which allow for extremely fast computation), our model regularly outperforms the alternative ones-which we reformulate here in terms of conditional models so as to allow for transparent comparison.

The paper is organized as follows. Sec. 2 reviews some preliminary notions concerning the Laplacian spectrum of graphs. The FRG model is presented in Sec. 3, where we also show how to estimate FRGs from data, and we analyze the dependence structures involved in FRGs. In Sec. 4 we review (and to some extent develop) a few alternative conditional approaches to random graph estimation, starting from some well-known generative models. All considered models are then evaluated experimentally in Sec. 5, while Sec. 6 draws some conclusions and sketches a few directions for further work.

GRAPHS, LAPLACIANS, AND EIGENVALUES

Let G = (V, E) be an undirected graph with n nodes.

In the following we assume that the graph is unweighted with adjacency matrix A. For each (unordered) pair of nodes {u, v} such that u = v, we take X uv to denote a binary random variable such that X uv = 1 if {u, v} ∈ E, and X uv = 0 otherwise. Since the graph is undirected, X uv = X vu .

The degree d u of a node u ∈ V is defined as the number of connections of u to other nodes, that is d u = |{v: {u, v} ∈ E}|. Accordingly, the degree matrix D of a graph G corresponds to the diagonal matrix with the vertex degrees d 1 , . . . , d n on the diagonal. The main tools exploited by the random graph model proposed here are the graph Laplacian matrices. Different graph Laplacians have been identified in the literature [START_REF] Von | A tutorial on spectral clustering[END_REF]. In this paper, we use consistently the unnormalized graph Laplacian, defined as L = D -A. The basic facts related to the unnormalized Laplacian matrix can be summarized as follows:

Proposition 1 [START_REF] Mohar | The Laplacian Spectrum of Graphs[END_REF]). In the following, the (algebraic) multiplicity of an eigenvalue λ will be denoted by M (λ, G). If the graph has one single connected component, then M (0, G) = 1, and the second smallest eigenvalue λ 2 (G) > 0 is called Fiedler eigenvalue. The Fiedler eigenvalue provides insight into several graph properties: when there is a nontrivial spectral gap, i.e. λ 2 (G) is clearly separated from 0, the graph has good expansion properties, stronger connectivity, and rapid convergence of random walks in the graph [START_REF] Mohar | The Laplacian Spectrum of Graphs[END_REF]. For example, it is known that λ 2 (G) ≤ µ(G), where µ(G) is the edge connectivity of the graph (i.e. the size of the smallest edge cut whose removal makes the graph disconnected). Notice that if the graph has more than one connected component, then λ 2 (G) will be also equal to zero, thus implying that the graph is not connected. Without loss of generality, we abuse the term Fiedler eigenvalue to denote the smallest eigenvalue different from zero, regardless of the number of connected components. In this paper, by Fiedler value we will mean the eigenvalue λ k+1 (G), where k = M (0, G).

For any pair of nodes u and v in a graph G = (V, E), we define two corresponding graphs G + and G -in the following way: G + = (V, E ∪ {{u, v}}), and G -= (V, E \ {{u, v}}). Clearly, we have that either

G + = G or G -= G.
A basic property concerning the Laplacian eigenvalues of G + and G -is the following [START_REF] Mohar | The Laplacian Spectrum of Graphs[END_REF][START_REF] Anderson | Eigenvalues of the Laplacian of a graph[END_REF], Cvetković et al., 1979]: v}}), and G -= (V, E \ {{u, v}}), then we have that:

Lemma 1. If G + and G -are two graphs with n nodes, such that {u, v} ⊆ V, G + = (V, E ∪ {{u,
1. n i=1 λ i (G + ) -λ i (G -) = 2; 2. λ i (G + ) ≤ λ i (G -) for any i such that 1 ≤ i ≤ n.

FIEDLER RANDOM GRAPHS

Fiedler random graphs are defined in Sec. 3.1, while in Secs. 3.2 and 3.3 we deal with the problems of estimating them from data and characterizing their statistical dependence structure respectively.

CONDITIONAL DISTRIBUTIONS

Given the notions reviewed above, we introduce the Fiedler delta function ∆λ 2 , which is defined as follows.

Let k = M (0, G + ). Then, ∆λ 2 (u, v, G) = λ k+1 (G + ) -λ k+1 (G -) (1)
In other words, for any pair of nodes u and v in graph G, the Fiedler delta value of the pair {u, v} in G is the (absolute) variation in the Fiedler eigenvalue of the graph Laplacian that would result from removing edge {u, v} from G + .

Lemma 1 immediately implies the following proposition:

Proposition 2. For any graph G = (V, E) and any pair of nodes u and v such that {u, v} ∈ E, we have that 0

≤ ∆λ 2 (u, v, G) ≤ 2. Proof. Let G -= (V, E \ {{u, v}}) and k = M (0, G).
The theorem follows straightforwardly from Lemma 1, given that ∆λ

2 (u, v, G) = λ k+1 (G) -λ k+1 (G -).
Using the Fiedler delta statistic, we can proceed to define a Fiedler random graph.

Given G = (V, E), let X (G) denote the set of random variables defined on G, that is X (G) = {X uv : u = v ∧ {u, v} ⊆ V}. Then: Definition 1. Given a graph G = (V, E), suppose we associate a subgraph G uv of G to each variable X uv ∈ X (G). We say that G is a Fiedler random graph if, for any X uv ∈ X (G), we have that P (X uv | X (G uv ) \ {X uv }) = P (X uv | ∆λ 2 (u, v, G uv )).
Concerning Definition 1, we stress the importance of two points. First, the key intuition motivating FRGs is to treat the Fiedler delta function as a real-valued random variable defined over graph configurations, and to exploit this random variable as a compact representation of those configurations. Second, FRGs are conditional models, i.e. they only define conditional distributions of edges in a random graph given the possible configurations of other edges. The model definition does not specify whether and how the considered conditional distributions can be combined to obtain a consistent factorization of the joint distribution of X (G).

In order to obtain that result, the model should take into account the global dependence structure of X (G), but this goal is not pursued by FRG modeling. Two interesting questions that can be raised on the nature of the FRG model are then the following. First, how does the Fiedler delta statistic behave with respect to the Markov dependence property [Besag, 1974, Frank and[START_REF] Frank | Markov Graphs[END_REF]? Second, how useful can be the information enclosed in the Fiedler delta statistic for learning and mining applications over large-scale networks? One basic result related to the first question is presented in Sec. 3.3, whereas Sec. 5 will address the second point.

MODEL ESTIMATION

The goal of estimating a FRG from data is to obtain an estimate of the conditional distribution P (X uv | X (G uv ) \ {X uv }), defined over subgraphs G uv of G. Starting from Definition 1, this quantity can be manipulated as follows:

P (X uv | X (G uv ) \ {X uv }) = P (X uv | ∆λ 2 (u, v, G uv )) = p(∆λ 2 (u, v, G uv )| X uv )P (X uv ) p(∆λ 2 (u, v, G uv )) = p(∆λ 2 (u, v, G uv )| X uv )P (X uv ) xuv p(∆λ 2 (u, v, G uv )| x uv )P (x uv ) (2)
In Eq. 2, P (X uv ) is the prior probability of observing an edge between any two nodes u and v, whereas

p(∆λ 2 (u, v, G uv )| X uv )
i , v i } in D. Let D xuv denote the set {u i , v i }: (x ui vi , G uivi ) ∈ D ∧ x uivi = x uv . Then, P (x uv ) = 1 n D xuv .
On the other hand, in order to estimate p(∆λ 2 (u, v, G uv )| X uv ), we need to make a choice concerning the form of the modeled density function. In the lack of prior knowledge concerning the form of the density, a reasonable choice is to adopt a nonparametric estimation technique. In particular, we use a kernel density estimator [START_REF] Rosenblatt | Remarks on Some Nonparametric Estimates of a Density Function[END_REF][START_REF] Parzen | On Estimation of a Probability Density Function and Mode[END_REF], defined as follows:

p(∆λ 2 (u, v, G uv )| x uv ) = = 1 |D xuv | {ui,vi}∈Dx uv K ∆λ2 {u, v}, {u i , v i } (3)
where

K ∆λ2 {u, v}, {u i , v i } = = 1 h K ∆λ 2 (u, v, G uv ) -∆λ 2 (u i , v i , G uivi ) h (4)
and K is a kernel function with bandwidth h [START_REF] Silverman | Density Estimation for Statistics and Data Analysis[END_REF]]. In our implementation of FRGs, we use the Epanechnikov kernel [START_REF] Epanechnikov | Nonparametric Estimation of a Multidimensional Probability Density[END_REF], namely

K(t) = 3
is more appealing than other possible functions (such as the Gaussian kernel) because of its computational convenience. At the same time, it displays analogous properties in terms of statistical consistency [START_REF] Silverman | Density Estimation for Statistics and Data Analysis[END_REF]].

DEPENDENCE STRUCTURE

In order to illustrate the behavior of FRGs with respect to conditional independence, we first recall the definition of Markov dependence for random graph models [START_REF] Frank | Markov Graphs[END_REF]

(X uv | X wz , N (X uv )) = P (X uv | N (X uv )).
Based on Definition 2, we say that the dependence structure of a random graph G is non-Markovian if, for disjoint pairs of nodes {u, v} and {w, z}, it does not imply that P

(X uv | X wz , N (X uv )) = P (X uv | N (X uv )), i.e. if it is consistent with the in- equality P (X uv | X wz , N (X uv )) = P (X uv | N (X uv )).
We can then prove the following proposition:

Proposition 3. There exist Fiedler random graphs with non-Markovian dependence structure.

Proof. To prove the proposition, we make use of the following equalities [START_REF] Fiedler | Algebraic connectivity of graphs[END_REF]: if graphs G 1 and G 2 are, respectively, a path and a circuit of size n, then

λ 2 (G 1 ) = 2 (1 -cos(π/n)) and λ 2 (G 2 ) = 2 (1 -cos(2π/n)).
Consider a graph G = (V, E) such that V = {u, v, w, z} and E = {{u, v}, {v, w}, {w, z}, {u, z}}, and let G uv = (V uv , E uv ) be the subgraph of G incident to E uv = E \{{w, z}}. If G is a Markov graph, it must be the case that P (X uv | X wz , X vw , X uz ) = P (X uv | X vw , X uz ). Now, suppose that G is a FRG. In this case, we have that

P (X uv | X wz , X vw , X uz ) = P (X uv | X (G) \ {X uv }) = p(∆λ 2 (u, v, G)| X uv )P (X uv ) p(∆λ 2 (u, v, G)) (5)
and, similarly,

P (X uv | X vw , X uz ) = p(∆λ 2 (u, v, G uv )| X uv )P (X uv ) p(∆λ 2 (u, v, G uv )) (6
) Eqs. 5-6 imply that, if the dependence structure of G is Markovian, then the following equality must hold:

p(∆λ 2 (u, v, G)| X uv ) p(∆λ 2 (u, v, G)) = p(∆λ 2 (u, v, G uv )| X uv ) p(∆λ 2 (u, v, G uv )) (7)
Since the configuration of G and G uv is given by a circuit and a path respectively, where both have size 4, we know that λ 2 (G) = 2 (1cos(π/2)) and

λ 2 (G uv ) = 2 (1 -cos(π/4)). Also, if G -= (V, E \ {{u, v}}) and G - uv = (V uv , E uv \ {{u, v}}), notice that λ 2 (G -) = λ 2 (G uv ), since G -is
also a path of size 4, and that M (0, G - uv ) = M (0, G uv ) + 1, since G - uv has one more connected component than G uv . Therefore, we have that ∆λ 2 (u, v, G) = 2 cos(π/4) and ∆λ 2 (u, v, G uv ) = 2 (1-cos(π/4)), i.e. ∆λ 2 (u, v, G) = ∆λ 2 (u, v, G uv ). Because of this inequality, there will exist parameterizations of p which do not satisfy Eq. 7, which means that the dependence structure of G is non-Markovian.

Note that the proof of Proposition 3 can be straightforwardly generalized to the dependence between X uv and X wz in circuits/paths of arbitrary size n, since the expression used for the Fiedler eigenvalues of such graphs holds for any n. In fact, suppose that the 4-nodes circuit G used in the proof is replaced with a circuit G * = (V * , E * ) of size n, where V * = V ∪ {s 1 , . . . , s 1 , . . . , s m , t 1 , . . . , t m } and E * is obtained from E by replacing {u, z} and {v, w}, respectively, with a path from u to z going through s 1 , . . . , s m and a path from v to w going through t 1 , . . . , t m , so that n = 2m + 4. In this case, if

G * uv is the subgraph of G * incident to E * uv = E * \ {{w, z}}, we have again that ∆λ 2 (u, v, G * ) = ∆λ 2 (u, v, G * uv ), which means that there exist FRGs such that p(∆λ2(u,v,G * )| Xuv ) p(∆λ2(u,v,G * )) = p(∆λ2(u,v,G * uv )| Xuv ) p(∆λ2(u,v,G * uv ))
.

OTHER CONDITIONAL RANDOM NETWORKS

In this section, we adapt three families of random graph models-namely the ERG, Watts-Strogatz (WS), and Barabási-Albert (BA) models respectively-to the task of conditional estimation of edge distributions, and we also review the entailed dependence structures. This will allow for a clear comparison of such models to FRGs in Sec. 5. Kronecker graphs are not considered here because it is not straightforward how they could be turned into conditional estimators.

EXPONENTIAL RANDOM GRAPH MODELS

Let G(V) be the set of all possible (undirected) graphs G with fixed vertex set V, where |G(V)| = √ 2 |V| 2 -|V| . Suppose that Φ denotes a collection of potential functions ϕ i over graphs, where each ϕ i (G) is a graph statistic such as the number of edges or the number of triangles in G. Then, an ERG with parameter vec-tor θ = (θ 1 , . . . , θ |Φ| ) defines the following Boltzmann distribution [START_REF] Garry | Recent developments in exponential random graph (p * ) models for social networks[END_REF]:

P (X (G)| θ) = 1 Z(θ) exp |Φ| i=1 θ i ϕ i (G) (8)
where each θ i is a real-valued parameter associated to the potential function ϕ i , and Z(θ) denotes the partition function [START_REF] Koller | Probabilistic Graphical Models: Principles and Techniques[END_REF], given by

Z(θ) = Gj ∈G(V) exp |Φ| i=1 θ i ϕ i (G j ) .
Computing the partition function exactly is clearly intractable for graphs with more than a few nodes, which makes it quite expensive to compute joint probabilities in large networks using ERGs.

ERGs entail a simple form for the conditional distribution of edges given the remainder of the graph:

P (X uv | X (G) \ {X uv }; θ) = exp |Φ| i=1 θ i ϕ i (G) xuv exp |Φ| i=1 θ i ϕ i (G xuv ) (9)
where G xuv denotes the configuration of G that we obtain by clamping the value of X uv to x uv . Therefore, if we want to estimate the parameters of an ERG so as to maximize the model log-likelihood with respect to a data sample D = (x u1v1 , G S1 ), . . . , (x unvn , G Sn ) , we can optimize each parameter θ i by taking the derivative ∂ ∂θi n j=1 log P (x uj vj | X (G Sj ) \ {X uv }; θ) and applying any gradient-based optimization technique.

Markov Random Graphs

Given a graph G = (V, E), define a k-star (for k ≥ 2) as a set E S k ⊆ E such that |E S k | = k and, for any {u, v} and {w, z} in E S k , {u, v} ∩ {w, z} = 1. Also, let a k-triangle (for k ≥ 1) be a set E T k ⊆ E of size 2k + 1 such that the elements of E T k are edges {u, v}, {u, w 1 }, {v, w 1 }, . . . , {u, w k }, {v, w k }, where w i = w j for any i and j such that 1 ≤ i, j ≤ k. If we use E(G), S k (G), and T k (G) to denote, respectively, the number of edges, the number of k-stars, and the number of k-triangles in G, then an exponential Markov random graph (MRG) model with parameter vector θ is defined by the following probability distribution [START_REF] Frank | Markov Graphs[END_REF]:

P (G| θ) = 1 Z exp η E(G) + K k=2 σ k S k (G) + τ T (G)
(10) where K is the maximum value that we want to consider for k in the k-star statistics, θ = (η, σ 2 , . . . , σ K , τ ), and T (G) = T 1 (G) is simply the number of triangles in G. That is to say, a MRG defines a Boltzmann distribution over graphs, with parameters corresponding to edge, k-star, and triangle statistics. Interestingly, the ER model can be characterized as a special case of MRG, where all the parameters except for η are set to zero.

MRGs are known to satisfy Definition 2 [START_REF] Garry | Recent developments in exponential random graph (p * ) models for social networks[END_REF]. Therefore, under the probability model defined by Eq. 10, we have that, for any pair of nodes {u, v} and any subgraph sample G S from G, P (X uv | X (G S ) \ {X uv }; θ) = P (X uv | N S (X uv ); θ), where N S (X uv ) is the set of all variables X wz from G S such that |{w, z}∩ {u, v}| = 1. Note that while Definition 2 specifies a general class of random graph models, MRGs in the strict sense refer to the class of ERGs defined above.

Higher-Order Models

Higher-order ERG models (HRGs) are readily obtained from MRGs by adding k-triangle counts (for k ≥ 1) to the Boltzmann distribution of Eq. 10 [Robins et al., 2007]. Moreover, in order to avoid fixing in advance the maximum value of the k parameter, a general formulation has been proposed for HRGs through the alternating k-star and the alternating k-triangle statistics [START_REF] Tom | New Specifications for Exponential Random Graph Models[END_REF], obtaining the following probability model:

P (G| θ) = 1 Z exp {η E(G) + σ S * (G) + τ T * (G)} (11)
where, if n is the number of nodes in G and S k and T k are the usual k-stars and k-triangle counts, S * (G) and T * (G) are defined as

S * (G) = n-1 k=2 (-1) k S k ρ k-2 and T * (G) = n-2 k=2 (-1) k T k
ρ k-2 (with ρ ≥ 1 acting as a sort of regularization parameter).

It is known that, under the distribution given by Eq. 11, P (X uv | X (G S ) \ {X uv }; θ) = P (X uv | N * S (X uv ); θ) [START_REF] Garry | Recent developments in exponential random graph (p * ) models for social networks[END_REF], where N * S (X uv ) is the set containing any X wz from G such that X wz = X uv and, for at least one edge {s, t}, we have that |{s, t} ∩ {u, v}| = 1 and {s, t} ∩ {w, z} = ∅. Clearly, N S (X uv ) ⊆ N * S (X uv ), which is why HRGs are 'higher-order' than MRGs.

WATTS-STROGATZ MODEL

The WS network model [START_REF] Watts | Collective dynamics of 'small-world' networks[END_REF]] defines a random network as a regular ring lattice which is randomly 'rewired' so as to introduce a certain amount of disorder, which typically leads to smallworld phenomena. Given nodes u 1 , . . . , u n , a WS network is generated by constructing a regular ring lattice such that each node is connected to exactly 2 δ other nodes. Network edges are then scanned sequentially, and each one of them is rewired with probability β, where, if i < j, rewiring an edge {u i , u j } means to replace it with another edge {u i , u k } such that k = i and u k is chosen uniformly at random from the set of all nodes that are not already linked to u i .

Interestingly, the degree distribution corresponding to a WS network G with parameters δ and β takes the following form [START_REF] Barrat | On the properties of smallworld network models[END_REF], for any degree k ≥ δ:

P (k) = I i=0 δ i (1 -β) i β δ-i (δ β) k-δ-i (k -δ -i)! exp(-β δ) (12
) where I = min{kδ, δ}. Given Eq. 12, we model the conditional distribution of a variable X uv given the remainder of G through the following quantity:

P (X uv | X (G) \ {X uv }) = P (d u (G)) P (d v (G)) xuv P (d u (G xuv )) P (d v (G xuv )) (13)
where d u (G) denotes the degree of node u in G, and each P is implicitly conditional on the values of δ and β. We refer to the conditional random graph model specified in Eq. 13 as the conditional Watts-Strogatz (CWS) model.

For the CWS model, the following proposition follows straightforwardly from Eq. 13:

Proposition 4. The dependence structure of CWS models is Markovian.

A simple estimate of the δ parameter for a data sample D = (x u1v1 , G u1v1 ), . . . , (x unvn , G unvn ) is the following:

δ = 1 2 n n i=1 d ui (G Si ) + d vi (G ui vi ) (14)
On the other hand, a simple strategy for estimating the rewiring probability by a maximum likelihood approach consists in parameterizing it as a sigmoid function β = 1 1+exp(-θ β ) , where θ β can be optimized by exploiting the derivative

∂ ∂θ β n i=1 log P (x uivi | X (G ui vi )\ {X uv }; δ, θ β ).

BARAB ÁSI-ALBERT MODEL

The BA model was originally proposed for explaining the scale-free degree distributions often observed in real-world networks [START_REF] Barabási | Emergence of scaling in random networks[END_REF]. In the BA model, the probability P (u) of linking to any particular node u in a network G = (V, E) takes the form Barabási, 2002], where α is a real-valued parameter affecting the shape of the degree distribution. Given P (u), we can use the following expression to characterize the conditional probability of observing edge {u, v} given the remainder of G [START_REF] Mark | Clustering and preferential attachment in growing networks[END_REF], Barabási et al., 2002]:

P (u) = du(G) α v∈V dv(G) α [Albert and
P (X uv = 1| X (G) \ {X uv }; α) = d u (G) α d v (G) α w∈V d w (G) α 2
(15) We refer to the random graph model of Eq. 15 as the conditional Barabási-Albert (CBA) model. The α parameter can be optimized straightforwardly by gradient-based methods, so as to maximize the objective n i=1 log P (x uivi | X (G uivi ) \ {X uv }; α) over an observed sample of size n.

The following property holds for CBA networks:

Proposition 5. The dependence structure of CBA models is non-Markovian.

Proof. Consider a CBA network G = (V, E), where {u, v} and {w, z} are edges in E such that {u, v} ∩ {w, z} = ∅. G is a Markov graph if and only if P (X uv | X wz , N (X uv )) = P (X uv | N (X uv )). Let G 1 and G 2 be the same as defined in the proof of Proposition 4. Although

d u (G 2 ) = d u (G 1 ) and d v (G 2 ) = d v (G 1 ), we have that w∈V d w (G 2 ) α > w∈V d w (G 1 ) α . Therefore, P (X uv | X wz , N (X uv )) = P (X uv | X (G 2 ) \ {X uv }) = d u (G 2 ) α d v (G 2 ) α w∈V d w (G 2 ) α 2 = d u (G 1 ) α d v (G 1 ) α w∈V d w (G 1 ) α 2 = P (X uv | X (G 1 ) \ {X uv }) = P (X uv | N (X uv )) (16)

EXPERIMENTAL EVALUATION

In order to assess experimentally the accuracy of FRGs as models of conditional edge distributions in largescale networks, we test them in the following link prediction setting. First, we take a network G = (V, E), with a number of nodes ranging from a few thousands to a couple of millions, and a number of edges going up to several millions. We then sample n training pairs {u i , v i } uniformly at random from V, with u i = v i .

For each sampled pair of nodes {u i , v i }, we recover the subgraph G uivi induced on G by the vertex set V uivi = {u i , v i }∪{w i : {u i , w i } ∈ E ∨{v i , w i } ∈ E}, so as to get a training set D = (x u1v1 , G u1v1 ), . . . , (x unvn , G unvn ) .

We then estimate from D each one of the conditional random graph models described in the previous sections, using the respective learning techniques.

Given the estimated models, we sample a test set T = (x u1v1 , G u1v1 ), . . . , (x umvm , G umvm ) from G such that T ∩ D = ∅, where each element of T is sampled through the same technique used for D. For each estimated model, we assess its accuracy by first using it to compute the conditional probability of observing a link for each pair of nodes {u j , v j } in T , and then calculating the area under the ROC curve (AUC) for the predictions delivered on T . Since the problem of predicting the presence of edges in real-world networks is severely unbalanced, i.e. |E| ≪ 1 2 | V| 2 -|V| , the ROC curve is a particularly appropriate evaluation measure [START_REF] Fawcett | An introduction to ROC analaysis[END_REF]. Also, in order to assess the sensitivity of each model to the size of the training sample, we ran some preliminary experiments with several different values of |D|. The result of this preliminary investigation was that, provided that the training sample contains at least a few pairs of linked nodes (in the order of ten or twenty), then the behavior of the different models is not significantly affected by increasing sample sizes. We compare the performance of the MRG, HRG, CWS, CBA, and FRG models on seven different networks, drawn from the arXiv e-print repository, the Enron email corpus, and the US road network respectively [START_REF] Leskovec | Graph evolution: Densification and shrinking diameters[END_REF][START_REF] Leskovec | Community Structure in Large Networks: Natural Cluster Sizes and the Absence of Large Well-Defined Clusters[END_REF]. The ER model is not considered in these experiments, since the involved independence assumption makes that model simply unusable for conditional estimation tasks. Some representative network statistics for the used datasets are summarized in than 0.5 (as recorded for some models) do not necessarily mean that prediction accuracy is being worse than random guessing. In fact, recent studies show that, depending on signal strength and stratification errors between training and test sets, random guessing can drop to AUC values even lower than 0.3 [START_REF] Parker | Stratification bias in low signal microarray studies[END_REF]. Therefore, AUC values less than 0.5 should be considered simply as no better than random guessing, rather than worse than random.

FRGs significantly outperform all other models both in the arXiv and in the Enron networks. Such networks are generally characterized by power law degree distributions, contrary to the road networks, which exhibit instead a regular topology, with degree distribution peaked around its mean. This explains the better behavior of the CBA model with respect to the remaining options in most of the scale-free networks. Interestingly, this trend is reverted for the HepPh data, where not only CWS outperforms CBA, but we also observe that HRG behaves much better than it does on the other scale-free networks. The fact that, in this particular case, FRGs remain the best available option provides evidence for their higher stability (i.e. robustness to variation in the underlying network properties) with respect to the other models. On the other hand, for the US road networks we observe a completely different pattern of results. This time, the FRG model performs worse than the CWS model or the two variants of ERG. The reason is related to the lattice-shaped regularity exhibited by such graphs. This is not surprising, as the WS model subsumes regular lattices as special case scenarios. A natural explanation for the superiority of the CWS model on the road networks lies in the fact that their parametric assumption concerning the observed degree distribution happens to be satisfied for such data, while it is clearly violated by the scale-free networks, which is exactly the opposite of what happens for the CBA model. It is worth noting, however, that FRGs are again performing significantly better than the CBA model. Therefore, FRGs appear to be a more stable class of conditional random network models, since their predictions ensure a regularly accurate behavior across networks displaying significantly different statistical properties. Overall, the presented results encourage the hypothesis that a genuinely nonparametric approach to conditional random graph modeling can offer dramatic advantages over parametric approaches.

CONCLUSION

The work presented in this paper started from two motivations. On the one hand, we remarked that statistical modeling of networks cries for nonparametric estimation, because of the inaccuracy often resulting from fallacious parametric assumptions. In this respect, we showed that statistics derived from the Laplacian spectrum of graphs (such as the Fiedler delta function) offer practical ways of developing nonparametric estimators. On the other hand, we suggested that a conditional approach to random graph modeling can be very effective in the large-scale setting, since it allows e.g. to dispense with the intractable partition functions often involved in joint distributions. With respect to this point, we showed how FRGs allow us to tackle very large datasets, processing effectively even millions of queries over networks with millions of nodes and edges.

Interesting options for future work consist in analyzing the effect of the subgraph sampling mechanism in determining the distribution of the Fiedler delta, as well as investigating alternative ways of exploiting the information enclosed in the graph spectrum.

Figure 1 :

 1 Figure 1: ROC curves for the AstroPh (a), GrQc (b), HepPh (c), and HepTh (d) networks.

Figure 2 :

 2 Figure 2: ROC curves for the Enron network.

Figure 3 :

 3 Figure 3: ROC curves for the RoadCA (a) and RoadTX (b) networks.

  . Let N (X uv ) denote the set {X wz : {w, z} ∈ E ∧ |{w, z} ∩ {u, v}| = 1}. A Markov random graph is then defined as follows:Definition 2. A random graph G is said to be a Markov graph (or to have a Markov dependence structure) if, for any pair of variables X uv and X wz in G such that {u, v} ∩ {w, z} = ∅, we have that P

Table 1 :

 1 General statistics for the network datasets used in the experiments. CC G and D G denote the average clustering coefficient and the network diameter respectively.

	Network AstroPh	|V| 18,772	|E| CC G 396,160 0.63	D G 14
	GrQc	5,242	28,980 0.52	17
	HepPh	12,008	237,010 0.61	13
	HepTh	9,877	51,971 0.47	17
	Enron	36,692	367,662 0.49	12
	RoadCA	1,965,206 5,533,214 0.04	850
	RoadTX	1,379,917 3,843,320 0.04 1,049

Table 1

 1

	. For the arXiv datasets

Table 2 :

 2 AUC values for the ROC curves plotted in Figs. 1-3. CC G and D G denote the average clustering coefficient and the network diameter respectively.

		AUC (%)
	Network CBA CWS FRG HRG MRG
	AstroPh 46.04 32.39 91.66 28.69 14.78
	GrQc	66.22 50.86 94.90 27.10 22.32
	HepPh	40.04 50.84 92.73 65.62 15.02
	HepTh	53.21 42.47 79.34 32.82 32.36
	Enron	48.45 49.03 85.92 19.94 11.23
	RoadCA 30.84 91.21 66.70 77.90 81.12
	RoadTX 34.67 93.01 63.75 62.51 68.36

1t 2 ½ 1 (|t|), where ½ y (x) = 1 if x ≤ y, and ½ y (x) = 0 otherwise. The Epanechnikov kernel

Acknowledgements

This work has been supported by the French National Research Agency (ANR-09-EMER-007). The authors are grateful to Romaric Gaudel, Rémi Gilleron, and Michal Valko for their suggestions and comments concerning this work.


Marc Bernard
email: marc.bernard@univ-st-etienne.fr

Elisa Fromont
email: elisa.fromont@univ-st-etienne.fr

Amaury Habrard
email: amaury.habrard@univ-st-etienne.fr

Marc Sebban
email: marc.sebban@univ-st-etienne.fr

Marc Sebban Handwritten Digit

Amaury Habard

Handwritten Digit Recognition using Edit Distance-Based KNN

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The Machine Learning (ML) course is given to the 1st year Master students at the University Jean Monnet in Saint Etienne. The course's duration is 60 hours for which 45 hours are dedicated to theory and desk exercises and 15 hours are dedicated to lab sessions on WEKA 1 or R 2 and on starting (∼ 8 hours) the student project on handwritten digit/letters recognition which also requires about 12 hours of personal work at home. These students already have a bachelor in computer science and thus have reasonable programming skills (for example they know how to use Java SWING 3 to create a nice GUI) and some experience in working 1 http://www.cs.waikato.ac.nz/ml/weka/ 2 http://www.r-project.org/ 3 http://docs.oracle.com/javase/tutorial/uiswing/ Teaching Machine Learning Workshop in small (4-8 persons) groups. They had a few hours of mathematics and in particular of statistics (∼ 60h) since high school but they cannot be considered as having a very strong mathematical and, in general, theoretical background.

The ML course covers topics related to supervised classification. The introductory part of this course (∼ 15h) covers the general definition of supervised classification, the error estimations, the non-parametric Bayesian methods and their link to the class estimations using Parzen windows and K-nearest neighbors (KNN). When focusing on the KNN part, a couple of hours are spent on the possible distances used to compare two instances and, in particular for stringbased data, on the edit distance. The other parts of the course cover different symbolic and numerical algorithms and models from Inductive Logic Programming to Hidden Markov Models through Decision Trees, Neural Networks etc. Given their background, this course is usually perceived as theoretical and the project is usually very well welcomed by the students.

The project is proposed to groups of up to 5 students. The goal of the project is to make students work on data acquisition and data storage, the implementation and optimization of a machine learning algorithm (here KNN and its cleaned and condensed versions), the evaluation of the parameters (e.g. the weights of the edit distance, the number of nearest neighbors, etc.) and of the system (cross validation, test set, etc.). A nice Graphical User Interface (GUI) is required to be able to give a live demo of the system at the end of the project.

Section 2 gives an overview of the project and of the theoretical prerequisites. Section 3 gives more details and a description of the steps the students go through when working on the project. We conclude in Sect. 4.

Project Description

In the following, everything which concerns a digit can also be applied to a letter.

Overview

The goal of the project is to design a handwritten digit recognition system based on a KNN classifier. At the end of the project, each group of students should be able to present a GUI on which a user can draw (for example with a mouse) a black and white digit on a screen. This drawing could be either, saved as an image (or directly as a Freeman code in a text file) to populate a database or, could be recognized online (and efficiently) by the system.

Freeman codes

There are multiple ways to encode a digit drawn on a screen for further processing. One of them consists in saving the digit and the frame in which it has been drawn as an image. From this image, a successful string representation has been proposed by Freeman [START_REF] Freeman | On the encoding of arbitrary geometric configurations[END_REF] to encode the digit. The computation of the Freeman codes is presented to the students during a lab session dedicated to the initiation of the project.

To compute the Freeman code of an image, one needs to move along a digital curve or a sequence of border pixels based on eight-connectivities. The direction of each movement is encoded by using a numbering scheme {i i = 0, 1, .., 7} denoting an angle of 45 * i (counter-clockwise) from the positive y-axis. The chain codes can be viewed as a connected sequence of straight-line segments with specified lengths and directions as showed in Figure 1.

KNN

We suggest the students to use the K-nearest neighbors algorithm for the digit recognition. From experience, we know that this is a simple yet effective method to automatically identify the drawn digit independently of the writing style or the scale (but the digit has to be drawn in a bounded-size window).

The KNN algorithm is presented during the lecture. Algorithm 1 explains how to classify an unknown example x using a dataset S. [START_REF] Vidal | New formulation and improvements of the nearest-neighbour approximating and eliminating search algorithm (aesa)[END_REF] for more details on the techniques).

• Let x be the example to classify by the nearestneighbor (NN) rule. Let us consider that the current NN of x is y which is at a distance δ from x.

• Let z be the next training example. If d(x, z) ≤ δ then update the current NN. Otherwise, remove the following two categories of examples located:

1. in the sphere centered at z and of radius d(x, z)δ, 2. out of the sphere centered at z and of radius d(x, z) + δ,

Distances

The distance used to compare two Freeman codes in the database is the edit distance (explained below).

A variant of the project could consist in testing the performance of other distances such as a L1, L2 or a Mahalanobis distance on other types of image representation (quad tree, 0-1 pixel matrix, interest point detectors,...).

Definition 1 The Levenshtein distance (or Edit Distance) between two strings x = x 1 ...x T and y = y 1 ...y V is given by the minimum number of edit operations needed to transform x into y, where an operation is an insertion, deletion, or substitution of a single character.

Rather than counting the number of edit operations, one can assign an edit cost to each edit operation and search for the less costly transformation:

• c(x i , y j ) is the cost of the substitution of x i by y j ,

• c(x i , λ) is the cost of the deletion (x i into the empty symbol λ),

• c(λ, y j ) is the cost of the insertion of y j .

Input:Two strings x(T ) and y(V ) Output: Edit Distance D(T, V ) between x and y D(0, 0) ← 0 for r=1 to T do

D(r, 0) ← D(r -1, 0) + c(xr, λ) end for for k=1 to V do D(0, k) ← D(0, k -1) + c(λ, y k ) end for for r=1 to T do for k=1 to V do d1 ← D(r -1, k -1) + c(xr, y k ) d2 ← D(r -1, k) + c(xr, λ) d3 ← D(r, k -1) + c(λ, y k ) D(r, k) ← min(d1, d2, d3) end for end for Return D(T, V )

Project Implementation

The first job of the group of students is to divide the different tasks between all the members of the group and organize the communication to be able to integrate all the parts before the evaluation of the project. A student version of the project can be found online at his address: http://labh-curien.univ-stetienne.fr/∼sebban/AA-RDF/applet.php?lang=en. The project focuses on machine learning. The aim is not to test the programming skills of the students. So, the choice of the programming language and the design of the GUI is left to the students. However, in our case, most of the students use the JAVA Swing API. This API provides tools to create a graphical user interface (GUI) for Java programs. An example of a resulting GUI is given Figure 2. The students have to create their own dataset and a structure to store the data. Some students decide to create a real database system (using for example Post-greSQL4 ) and store the images which contain a digit. With this kind of method, the student can use and compare multiple representations of the image. Some students decide to store only the Freeman code as textual information or as XML files as shown in Figure 3. Whatever the structure, the students find out early in their experimental process that more than one person should populate the dataset with some digits to make sure that the instances of each class are diverse enough (in term of writing style) to have a sufficiently good classification accuracy.

GUI

Dataset

Algorithms

The students can choose any programming language but it is usually done according to the choice of the API for the GUI (so most of the time, it is JAVA). Most of the necessarily algorithms are given (as pseudo code) during the lecture and have to be implemented. However, the part concerning the preprocessing of the digit (for example the transformation into a Freeman code independently of the position of the digit in the drawing window) is left as the students' choice. Depending on their organization skills, the students may not have time to implement all the possible optimizations of the KNN algorithm. In this case, we expect them to be critical about the performance of their algorithm especially during the demo.

Evaluation

Because of the limited size of their database, most of the students use a 5-fold or a leave-one-out cross validation strategy to evaluate the global performance of their algorithm. They keep increasing the size of their dataset until they reach a "sufficiently" good accuracy (around 80% with a representative enough database). Then, they try to use the dataset reduction techniques to evaluate the loss of accuracy.

They usually try to test different "K" for the KNN algorithm and keep the one that gives the best accuracy results by cross validation.

For the edit distance, they usually start by giving an equal weight to each of the three edit operations then they increase the weights of the insertion and deletion compared to the substitution. A better idea is to weight the substitution differently for each possible code. Indeed, it must seem easier to change a Freeman code 1 into a code 2 (because it is a small angle) than to change a Freeman code 1 into a code 4. [START_REF] Rico-Juan | Comparison of aesa and laesa search algorithms using string and tree-edit-distances[END_REF] is a good reference to find the weights of the edit operations for this problem.

Conclusion

The proposed project shows the whole process of supervised classification from the data acquisition to the design of a classification system and its evaluation. The project is usually extremely successful with the students. Most of them are very proud of the performance of their system and most of the time, they implement side games such as Sudoku or letter games (such as scrabble) to show the performance of their system online. The project, besides giving practical applications to some machine learning techniques seen in class, gives them the opportunity to work as a group. As explained in the text, some changes can be made on the representations of the images and on the distances used to compare two elements in the database.

Figure 1 .

 1 Figure 1. Representation of the digit 2 drawn on the left using a Freeman code.

Algorithm 1

 1 KNN algorithm Input:x, S, d Output: class of x for (x ′ , l ′ ) ∈ S do Compute the distance d(x ′ , x) end for Sort the |S| distances by increasing order Count the number of occurrences of each class l j among the k nearest neighbors Assign to x the most frequent class To be able to identify the digit online (with a real time constraint), the students have to optimize their original algorithm (e.g. by reducing the size of the database to reduce the computations). Two algorithms are proposed to speed up the classification process: one (see Algo. 2) which removes from the original dataset S the outliers and the examples located in the Bayesian error region, smoothing the decision boundaries, and one (see Algo. 3) which deletes the irrelevant examples (e.g. the centers of the homogeneous clusters). Algorithm 2 Step1 : Remove from the original dataset S the outliers and the examples of the Bayesian error region Input:S Output:S cleaned Split randomly S into two subsets S 1 and S 2 repeat Classify S 1 with S 2 using the 1-NN rule Remove from S 1 the misclassified instances Classify S 2 with the new set S 1 using the 1-NN rule Remove from S 2 the misclassified instances until stabilization of S 1 and S 2 S cleaned = S 1 ∪ S 2Another optimization can be made to speed up the original algorithm by speeding up the search for the nearest neighbor using the triangle inequality property

Figure 2 .

 2 Figure 2. Recognition of the digit 7 in the final GUI.

Figure 3 .

 3 Figure 3. Example of a dataset and its statistics computed by the students' system

http://www.postgresql.org/

Acknowledgements

This work is supported by the Pascal 2 network of Excellence.


Jérémy Robert
email: jeremy.robert@cran.uhp-nancy.fr

J.-P Georges

Éric Rondeau

Thierry Divoux

Jean-Philippe Georges

Minimum Cycle Time Analysis of Ethernet-Based Real-Time Protocols

Keywords: Real-time systems, performance analysis, embedded systems

The Ethernet standard is a standard solution for interconnecting industrial devices despite its intrinsic drawbacks, particularly its nondeterministic medium access method. Many Ethernet-based commercial solutions available (COTS -Components Off the Shelves) on the market guarantee time performance. This means that user selection of one particular solution is a critical decision, but the choice often depends more on political strategizing with an industrial device manufacturer than on the intrinsic performance of Ethernet-based interfaces. The objective of this paper is to provide a formal behavioural analysis of each Ethernet-based solution, in order to facilitate comparison. We also consider the issue of synchronization of distributed nodes.

Introduction

Fieldbuses interconnecting industrial equipment are typically designed by Programmable Logic Controller (PLC) manufacturers. This is mainly attributable to specific constraints on industrial communication, which require a high degree of expertise. Fieldbuses must be robust in the noisy environments produced by the plant (physical layer), deterministic in guaranteeing data refresh during controller cycle periods (link layer), and capable of exchanging information between all types of industrial devices (application layer). Manufacturers provide different solutions to satisfy the first two constraints, without necessarily considering the final constraint to be relevant in business terms. This has resulted in the specification of a large number of fieldbus standards as described in [START_REF] Felser | The fieldbus war: History or short break between battles?[END_REF].

A new trend endorsed particularly by the IAONA (Industrial Automation Open Networking Alliance) consortium was to promote the Ethernet network as a standard for industrial communications. The expected benefits are less costly network installations, because equipment is available off-the-shelf, and the avoidance of interoperability problems, because Ethernet technology is broadly used. Other advantages are that Ethernet is a well-known protocol, which is widely implemented, and its performance improves continuously with technological evolution (especially bandwidth).

However, access to the Ethernet medium relies on the nondeterministic CSMA/CD algorithm, which applies a stochastic method to resolving collisions and cannot guarantee that message transmissions will be received in bounded time. Consequently, the native Ethernet protocol cannot be directly implemented in a plant with severe time constraints, and many Ethernet-based solutions have been proposed to overcome this issue, as mentioned in [START_REF] Alves | Ethernet goes real-time: a survey on research and technological developments[END_REF][START_REF] Decotignie | A perspective on ethernet-tcp/ip as a fieldbus[END_REF][START_REF] Decotignie | Ethernet-based real-time and industrial communications[END_REF]. However, if these Ethernet-based solutions are adapted to industry requirements, they lead to two types of problem. The first is that different solutions are not interoperable because different fieldbuses are developed by each manufacturer. The second problem is that time performances are insufficiently evaluated or compared.

According to [START_REF] Jasperneite | Limits of increasing the performance of industrial ethernet protocols[END_REF], different Ethernet products can be summarily classified into three main categories: the native Ethernet standard (Ethernet/IP, Modbus/TCP), Ethernet solutions using the priorities defined in IEEE802.1D/Q, and Ethernet-based solutions that incorporate new scheduling features in ASIC/FPGA (EtherCAT, Profinet IRT).

The last approach enables the elimination of all collisions and simplifies transmission time estimation, as described in [START_REF] Jasperneite | Limits of increasing the performance of industrial ethernet protocols[END_REF][START_REF] Prytz | A performance analysis of ethercat and profinet irt[END_REF]. The authors compared EtherCAT with Profinet IRT in a simplified context using analytic models. This analysis was improved and refined in section 2 and extended to two other well-known industrial Ethernet products: the Modbus/TCP solution and Ethernet/IP. In section 3, we present several scenarios that use the analytic models to facilitate assessment of different Ethernet products.

2

Estimation of minimum cycle time

Introduction

The general objective of this study was to compare the time performance of the major industrial Ethernet products available on the market. This comparison could be achieved only in a common application context. Thus, the specification of the communication scenario was based on one controller (for example, a PLC), interconnecting sensors, and actuators in an Ethernet network.

The controller was treated as the communication master that initiated all dialogues with slave nodes (sensors and actuators). The controller was characterized by its controller cycle time period, which was divided into three steps, as shown in Figure 1: -the sensor data refresh time, in the controller memory, -the processing time, and -the actuator update time.

Steps (1) and (3) represent communication periods, which should be less than the difference between the controller cycle time period and the processing time period. Thus, the time performance of each Ethernet product was compared according to a constraint named minimum cycle time, which was defined as:

The minimum cycle time was the communication time required by the controller to both collect and update the data memories of all sensors and actuators. In the following section, the analytic models of minimum cycle time are elaborated for EtherCAT, Profinet IRT, Modbus/TCP, and Ethernet/IP. These models all used the following parameters: the transmission delay, the network device latency, the propagation delay, the link capacity, the payload, and the number of slaves. The notations of these are given in table 1. It was assumed that there were no transmission errors and that the network was dedicated to the PLC application and was not shared with other applications.

EtherCAT

The EtherCAT network was developed by the Beckhoff company (type 12 in standard IEC 61158, [START_REF] Iec | Digital data communications for measurement and control ? fieldbus for use in industrial control systems: Part 3: Data link service definition[END_REF][START_REF] Iec | Digital data communications for measurement and control ? fieldbus for use in industrial control systems: Part 4: Data link protocol specification[END_REF]). In theory, EtherCAT cards are standard Ethernet interfaces. In practice, specific hardware (FPGA, Field-Programmable Gate Array, or ASIC, Application-Specific Integrated Circuit) is used to mitigate the frame forwarding delay. The EtherCAT network adds a master/slave protocol over the Ethernet. A frame is sent by the master and slaves can read and write data on the fly. The duration of reading or writing operations corresponds only to the network device latency ℓ, which is independent of frame size and the same for all slaves. A logical ring is defined between the slaves such that when a frame reaches the last slave in the ring it is returned to the master via all the slaves. The space-time diagram shown in Figure 2 illustrates the behaviour of EtherCAT communications. The EtherCAT protocol can support both line and ring topologies. Because the line topology is used mainly in the industrial framework, this topology was investigated in our study (Figure 3). The EtherCAT datagram is directly encapsulated inside the basic Ethernet frame as shown in Figure 4. An EtherCAT frame is composed of a header specifying the length of the frame and a list of datagrams. The number of datagrams depends on the number of slaves. A datagram is defined for each slave, and it contains the command type and associated data. In this study, our analysis of EtherCAT performance considered the following hypothetical test scenario: -the topology was a line, -the initialization step was ignored and only cyclic communication was studied, -the master sent only one frame per cycle, and -the payload x was the same for each slave.

The minimum cycle time shown in Figure 2 was determined for this scenario.

The link transmission delay is the ratio between frame size and link capacity C. The total frame size can be divided into two parts: -a constant value that equals the sum of the Ethernet protocol (26 bytes), the interframe gap (corresponding to the time of 12 bytes), and the Ether-CAT header (2 bytes) and -a variable value that depends on the slave number n, the amount of data to transport x, and the header (12 bytes).

The link transmission delay is:

τ = 8 (40 + max (44, n (12 + x))) C . (1)
The term 44 in the equation ( 1) was added to ensure the minimum data size defined by the Ethernet protocol. If the EtherCAT telegram length was less than 46 bytes, an equivalent amount of padding was inserted in the Ethernet frame. The EtherCAT telegram already included a 2 bytes header, which meant that there was no padding requirement when the length of the datagram sequence was larger than 44 bytes. As shown in Figure 2, the cycle time was estimated using the expression:

Γ = (2n -1) ℓ + 2nδ + τ = (2n -1) ℓ + 2nδ + 8 (40 + max (44, n (12 + x))) C . (2)
It should be noted that (2) considers only one frame. Because the Ethernet payload size depends directly on the number of slaves and the Ethernet frame size cannot exceed 1 526 bytes (and therefore the data field 1 500 bytes), (2) is only valid if the number of devices interconnected to the network is less than:

n ≤ n max = 1500 -EH 12 + x ,
where EH is the EtherCAT header size (2 bytes) and n max is the maximum number of datagrams, of length x, which can be included in a single frame (in the following, we have assumed that x ≤ 1486 bytes).

In general, the number of slave devices on the network can be greater than the frame size capacity. This means that the controller has to send more than one frame in a cycle time. In fact, the number of Ethernet frames required to support n devices with a constant payload x is given by:

k = n n max .
Consequently, (2) now integrates a different transmission time for each frame and finally gives:

Γ = (2n -1) ℓ + 2nδ + 8 C 40k + (k -1) n max (12 + x) + 8 C max 44, (n -(k -1) n max ) (12 + x) (3)
The final term of expression (3) was used to differentiate cases where the last frame generated padding.

Similar results are given in [START_REF] Jasperneite | Limits of increasing the performance of industrial ethernet protocols[END_REF] but, in contrast to this earlier work, expression (3):

• takes into account the on the fly minimum cycle time mechanism proposed by EtherCAT; the main advantage of this is that a device can begin frame forwarding before complete reception of a frame (in contrast to, e.g., store-and-forward mode), which significantly reduces the forwarding time as shown in Figure 2,

• considers the use of padding, as defined by Ethernet,

• integrates the time required to forward the information sent from devices to the controller, as shown in Figure 2, and

• considers cases where the number of slaves and their payload requires the utilization of several frames.

The accuracy of the EtherCAT synchronization mechanism was reported in [START_REF] Cena | On the accuracy of the distributed clock mechanism in ethercat[END_REF], which shows that this issue need not be considered because it was estimated as equal to a few nanoseconds.

Profinet IRT

Introduction

The Profinet protocol was developed by the Siemens company (type 10 in standard IEC 61158, [START_REF] Iec | Digital data communications for measurement and control ? fieldbus for use in industrial control systems: Part 3: Data link service definition[END_REF][START_REF] Iec | Digital data communications for measurement and control ? fieldbus for use in industrial control systems: Part 4: Data link protocol specification[END_REF]). Profinet IRT manages real-time communications. However, standard Ethernet cards cannot be used because Profinet IRT requires the operation of specific hardware on slaves (ASIC type -2 or 4 ports inbuilt switch). Profinet IRT is based on the time-slice mechanism, which specifies two modes, the asynchronous mode and the isochronous mode

The asynchronous and isochronous modes are used for unconstrained traffic and real-time traffic, respectively. Our study dealt only with real-time traffic. Thus, only the isochronous mode was analysed. For more information on these two modes, the reader can refer to [START_REF] Jasperneite | Investigations on a distributed time-triggered ethernet realtime protocol used by profinet[END_REF]. In the isochronous mode, the master sends one data frame to each device and each device replies to the master. Profinet can support line, star, and ring topologies. Only the line topology (see Figure 5), with full-duplex links, was analysed in this study (as for EtherCAT).

In this study, the Profinet IRT used the slipstreaming effect, where the controller began by sending frames to the most remote slave in the line topology, and then to the second remote slave slave and so on, until it reached the nearest slave. This mechanism enabled a reduction in the cycle time by minimizing the transmission time.

The slipstreaming effect was also applied to exchanges from slaves to the controller. The links were set up in full-duplex mode. The global communication scheme is given in Figure 6. The space-time diagram shown in Figure 6 illustrates the optimum use of the Profinet IRT protocol. Devices were assumed to be synchronized using a clock synchronization protocol, such as the IEEE 1588 standard. The IEEE 1588 standard generates synchronization frames (PTP) between devices, but these frames were not relevant to this study.

Minimum cycle time estimation

The slave devices periodically sent their messages to the controller at the same times as messages were sent by the controller to slaves.

Only the positive characteristics of the slipstream effect were considered. This required that τ ≥ δ + ℓ, as noted by [START_REF] Jasperneite | Limits of increasing the performance of industrial ethernet protocols[END_REF].

In such cases, the minimum cycle time was given by [START_REF] Jasperneite | Limits of increasing the performance of industrial ethernet protocols[END_REF][START_REF] Prytz | A performance analysis of ethercat and profinet irt[END_REF] as the sum of:

• the latency ℓ that crosses all devices plus the propagation delay δ for each link and

• and the link transmission time τ for each frame sent by the controller.

Hence, the minimum cycle time is written as:

Γ = δ + ℓ + nτ. (4)
Equation ( 4) was developed by analysing the transmission delay τ . The Profinet datagram is encapsulated in the Ethernet frame as shown in Figure 7.

Pre. Three fields are added to the Ethernet frame: a data identifier (2 bytes), the data value (with the assumption that x ≤ 1494 bytes), and an information status (4 bytes). When considering the constraint of the minimal Ethernet frame size, the transmission delay of a Profinet frame was given by:

τ = 8 38 + max (46, 6 + x) C . (5)
The term 38 in the equation ( 5) corresponds to the size of the Ethernet layer (26 bytes) plus the interframe gap (12 bytes). As with EtherCAT, the term 46 was added in order to ensure the minimum data size as defined by the Ethernet protocol. If the EtherCAT telegram length was less than 46 bytes, an equivalent amount of padding was inserted in the Ethernet frame.

The final form of equation ( 4) was expressed as:

Γ = δ + ℓ + n 8 C 38 + max (46, 6 + x)

Comments

The EtherCAT minimum cycle (3) was less than the Profinet IRT one (4), because the Profinet IRT transmission delay was multiplied by the number of devices. The Profinet IRT analysis assumed that all devices had the same clock reference. Figure 6 shows that all devices were synchronized because of the IEEE 1588 protocol. They shared the same clock and periodically sent their messages at the same time. Hence, it may be expected that clock synchronization errors will increase the minimum cycle time. This paper aims only at comparing optimal performances of COTS Ethernet-based protocols, i.e. Profinet IRT nodes sharing the same clock reference. Next studies will hence aim at reporting this synchronization issue.

Modbus/TCP

Modbus is a serial communication protocol developed by Modicon in 1979. Modbus/TCP is a variant of the Modbus protocol (type 15 in standard IEC 61158, [START_REF] Iec | Digital data communications for measurement and control ? fieldbus for use in industrial control systems: Part 3: Data link service definition[END_REF][START_REF] Iec | Digital data communications for measurement and control ? fieldbus for use in industrial control systems: Part 4: Data link protocol specification[END_REF]), which uses the Ethernet physical and link layers [START_REF]Modbus-IDA[END_REF]. Modbus/TCP encapsulates a Modbus frame into a TCP frame as shown in Figure 8. The Modbus datagram is composed of the ModBus Application Protocol Header (MBAP), the function (read/write), and the data value.

Pre. Modbus/TCP is a pragmatic approach that works on several types of configurations that impact on its performance. Parameters, including topology, exchange management in the application layer, and the processing capacity of devices, change the Modbus/TCP time behaviour.

Modbus/TCP is based on connection-oriented transactions and it can use different exchange models, such as master/slaves, producer/consumer, or client/server. The objective of this study was to compare several Ethernet products, so it was necessary to use similar contexts for all Ethernet products. Because EtherCat and Profinet IRT were previously analysed using a master/slaves model, this model was also used for Modbus/TCP evaluation. The master/slaves model also simplified the analysis because all the devices were synchronized by network events (none clock synchronization protocol is required).

The communication scheme defined in this study followed these steps. The master sent a frame to one slave and when the slave received this frame, it sent a reply to the master. When the master received the reply, it repeated the same procedure with another slave. All the slaves were processed by the master using a round-robin method. This communication scheme was implemented in the application layer of the OSI model. The request/reply protocol is shown in Figure 9. Modbus/TCP not only supports application data but also TCP PDU, and it comprises opening and closing TCP connections when acknowledging segment reception. The acknowledgment can be achieved either immediately a segment is received, after the reception of several segments, or inside the next data transmission (piggybacking). In practice, TCP behaviour changes according to the operating system, the TCP configuration, and whether or not the Nagle algorithm is used. In this study, we assumed that the timeout to send the acknowledgement was 0.5 s, meaning that only piggybacking was analysed, as shown in Figure 9. The transient states of the TCP opening and closing steps were not considered in the modelling.

Modbus/TCP supports line, star, and ring topologies. Complex architectures based on switches can be used, especially when the network is shared by several applications, although switches induce additional costs. However, because only the master/slaves communication scheme was considered in this study, a bus infrastructure was used It was assumed that each device was interconnected through a common hub, and the propagation delay was equal to δ. Figure 9 shows that the cycle time was equal to the number of devices n multiplied by the time required to poll one device.

Given previous assumptions, this can be written as:

Γ = n (τ req + τ rep + 2 (2δ + ℓ)) . (6)
In the worst case, where a TCP acknowledgment is sent for each segment, the minimal cycle time would be equal to:Gamma = n (τ req + τ rep + 3 × (2δ + ℓ) + 2τ ack ) where τ ack = 672/C and 672 is the minimal size in bits for an Ethernet frame (interframe gap included). When the architecture was composed of several hubs, the propagation delay increased and ( 6) was slightly different.

The frame transmission time was composed of:

• a constant part related to the sum of the Ethernet protocol (38 bytes with the interframe gap), the IP header (20 bytes without options), the TCP header (20 bytes without options), and the Modbus/TCP header (7 bytes for the ModBus Application Protocol Header),

• a variable part related to the type of Modbus message, with the size changing according to the type of data (function code) and the transaction model state (request or reply).Consequently, the analysis of Modbus/TCP was only valid for an application framework. In this study, we considered only write requests.

• a variable part proportional to the payload (because the byte count field was stored as a single byte, this indicates that the date field length was limited to x ≤ 255 bytes).

The sum of IP, TCP, and Modbus header sizes was larger than the minimal Ethernet data length. Thus, no added padding was required. The delay is directly given by:

τ req10h + τ rep10h = 8 91 + x C + 8 90 C
such that (6) corresponds to:

Γ 10h = n 8 181 + x C + 2 (2δ + ℓ) . (7)
2.5 EtherNet/IP

Introduction

EtherNet/IP (IP, Industrial Protocol) is a network developed by Rockwell Automation in 2001 and supported by ODVA (Open DeviceNet Vendor Association) [START_REF] Brooks | Ethernet/IP-industrial protocol[END_REF][START_REF] Schiffer | The common industrial protocol (CIP) and the family of CIP networks[END_REF]. EtherNet/IP (type 2 in standard IEC 61158, [START_REF] Iec | Digital data communications for measurement and control ? fieldbus for use in industrial control systems: Part 3: Data link service definition[END_REF][START_REF] Iec | Digital data communications for measurement and control ? fieldbus for use in industrial control systems: Part 4: Data link protocol specification[END_REF]) uses the Common Industrial Protocol (CIP) for off-the-shelf Ethernet products and TCP-UDP/IP stack. Ethernet/IP is a connection-based network. A CIP connection defines the type of packet sent to the network. There are two types of connections: the Explicit Messaging connection and the I/O (or Implicit) connection. Explicit Messaging provides generic and multi-purpose communication paths between two nodes, whereas I/O messaging is specific to application I/O data and provides serial purpose communication paths. When the application is time-constrained, I/O Messaging is the preferred mode because it employs UDP rather than TCP sockets. CIP uses the producer/consumer model and requires broadcast exchanges encapsulated in UDP. This study evaluated only I/O connections.

Because EtherNet/IP relies on COTS, no particular topology was specified. Bus-or switchbased architectures are both possible. Switches are interesting because they break the collision domain, allowing the support of VLAN and the classification of service mechanisms [START_REF] Modlovansky | Utilization of modern switching technology in ethernet/IP networks[END_REF]. Switched architectures are recommended for exchange management with time-critical (implicit) messaging. However, several switched Ethernet architectures are possible. In this study, a linear switched topology was selected, as shown in Figure 10. This can be viewed as an extension of the experimental set-up considered in [START_REF] Alessandria | Performance analysis of ethernet/ip networks[END_REF], where the number of switches varied according to the number of ports per switch.

In contrast to Modbus/TCP, utilization of a slipstreaming effect on Ethernet/IP does not facilitate the use of a bus topology. Indeed, it is possible that two messages can be on the network at the same time (as shown in Figure 11).

In order to compare the other protocols with Ethernet/IP, a similar exchange scenario was proposed. The controller sent a frame to each slave device and each slave device produced data that were sent to the controller. Figure 11 Initially, EtherNet/IP did not support medium access synchronization points, as found in the master/slaves technique used in Modbus/TCP. This meant that any device could access the network at any time. However, the time synchronization mechanism can be used because of the support of EtherNet/IP by the IEEE 1588 [START_REF]IEEE standard for a precision clock synchronization protocol for networked measurement and control systems, ANSI/IEEE standard[END_REF] protocol (implemented by the CIPSync profile). Consequently, devices can send messages using the same clock reference. We used the slipstreaming effect for the controller and a common departure time for the devices, as shown in Figure 11. Obviously, this profile corresponded to an ideal case related to the minimal cycle time for this architecture, because lock synchronization errors would lead to another profile inducing a longer cycle time and the use of the slipstreaming effect requires (see Profinet IRT) that τ ≥ δ + ℓ.

Minimum cycle time

Considering the profile given in Figure 11, the minimum cycle time corresponded to the sum of: -the latency ℓ for crossing only one switch plus twice the propagation delay δ, between a controller/device and a switch and -the link transmission time τ for each frame sent by the controller.

Thus, the minimum cycle time can be written as:

Γ = 2δ + ℓ + nτ. (8)
The term given in [START_REF] Iec | Digital data communications for measurement and control ? fieldbus for use in industrial control systems: Part 4: Data link protocol specification[END_REF] was valid only if τ ≥ δ + ℓ, meaning that switches did not work in the store-and-forward mode. The link transmission time was computed according to the encapsulation format of I/O messages Figure 12 shows that 18 bytes were added by the CIP protocol, 8 bytes by UDP, 20 bytes by IP, and finally 38 bytes by Ethernet. In this case, there was no padding. Thus, under the (8) can be modified as:

Γ = 2δ + ℓ + 8n 84 + x C . (9)
The profile given in Figure 11 can be optimized if several "items" (data) are encapsulated inside one frame, as specified in the CIP protocol. It induces a mitigation of the minimum cycle time. Since we employed the same approach that for Profinet IRT (in particular the slipstreaming effect), next studies will hence aim at reporting the synchronization effect on the minimum cycle time.

the next section, the results obtained for EtherCAT, Profinet IRT, Modbus/TCP, and EtherNet/IP are compared in different application contexts.

Comparisons

The objective was to analyse and compare the behaviour of different Ethernet-based solutions. The topologies used in each assessment were linear. Two bandwidths were studied: 100 Mb/s and 1 Gb/s. We analysed two payload sizes introduced by [START_REF] Prytz | A performance analysis of ethercat and profinet irt[END_REF] as being representative of an industrial context: 16 bytes and 100 bytes. The network device latencies ℓ were defined in [START_REF] Jasperneite | Limits of increasing the performance of industrial ethernet protocols[END_REF][START_REF] Prytz | A performance analysis of ethercat and profinet irt[END_REF] and those used in our study are described in Table 2. It was assumed the performances of switches used in EtherNet/IP were the same as those used by Profinet IRT. Clearly, alternative assumed values would yield different results. Thus, the results for the study using EtherNet/IP are only valid when τ ≥ δ + ℓ. The link propagation delay was 50 ns and corresponded to a distance of 10 m between two devices. The comparison is related to the minimum cycle time without synchronization errors. The first observation was that Modbus/TCP provided the worst results, whatever the payload. This was because of its medium access mechanism, which is based on polling at the application level.

In the case of a small payload, EtherCat provided the best results in the FastEthernet mode. The impact of bandwidth on Profinet IRT and Ethernet/IP was very significant, because it enabled a reduction in the minimum cycle time. This benefit was lower with EtherCat. The explanation is simple: EtherCat sends only one frame for communicating with all its slaves and so EtherCat provides only one link per transmission. In contrast, Profinet IRT and Ethernet/IP sent n frames to dialogue with n slaves, which provides n times the link transmission. The bandwidth was the most crucial parameter in reducing the time cycle. Moreover, 2n latencies must be considered with EtherCat because Profinet and EtherNet/IP cycle times were composed of only one switch latency ℓ, because of the effect. Profinet IRT provided the best results at 1 Gb/s.

Figure 13 shows that all the Ethernet solutions, with the exception of FastEthernet Modbus/TCP, could interconnect with more than 60 slaves in cycle times of less than 1 ms. This cycle time constraint, and the number of slaves and the frame size, enables adequate coverage for most industrial applications. Performances were similar when the payload was increased. The only difference was that the EtherCat minimal cycle time was larger than Profinet when the number of devices was increased. This was because the EtherCat telegram had to be fragmented, which decreased protocol performance. An increase in frame size means that the choice of the FastEthernet solutions must be carefully considered because the cycle time grows quickly. The GigaEthernet solution eliminated this problem.

Conclusion

This study analysed the time performance of Industrial Ethernet protocols. The general conclusion was that all Ethernet protocols perform suitably with a bandwidth of 1 Gb/s for in-terconnected real-time systems. At 100 Mb/s, special attention is required by engineers in the selection of Ethernet-based protocols. Other considerations have to be taken into account in the selection of Ethernet-based products, including persistence in the market and interoperability with other industrial equipment and network. Regarding interoperability concerns, Ethernet/IP is a standardized solution and provides acceptable performance for controlling industrial systems. Ethernet/IP can also implement priority mechanisms (IEEE 802.1p), which is important when a network is shared with other applications because this facility enables the differentiation of network services offered by real-time traffic and unconstrained traffic.

Future works are related to the clock synchronization issue, especially for solutions based on the slipstreaming effect. The objective is to deal with non optimal situations for Profinet IRT and Ethernet/IP where synchronization errors appear.

Figure 1 :

 1 Figure 1: PLC cycle time

Figure 2 :

 2 Figure 2: EtherCAT space-time diagram

Figure 3 :

 3 Figure 3: EtherCAT line topology

Pre. 8 DA 6 SA 6 ET 2 EH 2 DG 1 .

 866221 , SA/DA: Source/Destination Address, ET: Ethertype EH: EtherCAT Header, DG: Datagram, H: Header, Wc: Working counter

Figure 4 :

 4 Figure 4: EtherCAT frame (field lengths are given in bytes)

  Figure 5: Profinet line topology

Figure 6 :

 6 Figure 6: Profinet IRT space-time diagram with slipstreaming effect

Figure 7 :

 7 Figure 7: Profinet IRT frame (field lengths are given in bytes)

2 Figure 8 :

 28 Figure 8: General Modbus/TCP frame and function specific Modbus application frame (field lengths are given in bytes)

  controller

Figure 9 :

 9 Figure 9: Modbus/TCP space-time diagram

Figure 11 :

 11 Figure 11: EtherNet/IP space-time diagram for 4-ports switches

2 T 2 AL 2 Conn. 4 4 DT 2 DL 2 dataxFigure 12 :

 222442212 Figure 12: EtherNet/IP I/O messaging frame (field lengths are given in bytes)

Figure 13 :

 13 Figure 13: Minimum cycle time as a function of the number of devices with a constant payload of 100 bytes per device and two bit rates of 100 Mb/s and 1 Gb/s

Table 1 :

 1 Notations

	Terms	Notation Units
	Minimum cycle time	Γ	s
	Delivery time (from frame i)	d i	s
	Transmission delay	τ	s
	Network device latency	ℓ	s
	Propagation delay	δ	s
	Link capacity	C	bits/s
	Payload	x	bytes
	Number of network devices (slaves)	n	-

  shows this behaviour.

	device 1	device 2	device n
	controller
	Figure 10: EtherNet/IP switched linear topology
	controller	t
	switch A	t
	device 1	t
	device 2	t
	switch B	t
	device 3	t
	communication cycle time

Table 2

 2

		: Latencies
	Protocol	FastEthernet (100 M b/s) GigaEthernet (1 Gb/s)
	EtherCat	1.35 µs	0.85 µs
	Profinet IRT	3 µs	0.6 µs
	Modbus/TCP	1 µs (hub)
	EtherNet/IP	3 µs	0.6 µs
	13 show the minimum cycle times (in ms), according to the number of slave devices with a
	payload equal to 100 bytes