Konstantin Avrachenkov
email: k.avrachenkov@sophia.inria.fr

Natalia M Markovich
email: markovic@ipu.rssi.ru

Jithin K Sreedharan
email: jithin.sreedharan@inria.fr

Distribution and Dependence of Extremes in Network Sampling Processes

Keywords: Network sampling, extreme value theory, extremal index, random walks on graph

We explore the dependence structure in the sampled sequence of large networks. We consider randomized algorithms to sample the nodes and study extremal properties in any associated stationary sequence of characteristics of interest like node degrees, number of followers or income of the nodes in Online Social Networks etc, which satisfy two mixing conditions. Several useful extremes of the sampled sequence like kth largest value, clusters of exceedances over a threshold, first hitting time of a large value etc are investigated. We abstract the dependence and the statistics of extremes into a single parameter that appears in Extreme Value Theory, called extremal index (EI). In this work, we derive this parameter analytically and also estimate it empirically. We propose the use of EI as a parameter to compare different sampling procedures. As a specific example, degree correlations between neighboring nodes are studied in detail with three prominent random walks as sampling techniques.

Mots-clés : Échantillonnage de réseau, théorie de valeur extrême, indice extrémal, marche aléatoire sur un graphe.

Introduction

Data from real complex networks shows that correlations exist in various forms, for instance the existence of social relationships and interests in social networks. Degree correlations between neighbors, correlations in income, followers of users and number of likes of specific pages in social networks are some examples, to name a few. These kind of correlations have several implications in network structure, for example, degree-degree correlations manifests itself in assortativity or disassortativity of the network [START_REF] Barrat | Dynamical processes on Complex Networks[END_REF].

We consider very large complex networks where it is impractical to have a complete picture a priori. Crawling or sampling techniques are in practice to explore such networks by making use of API calls or HTML scrapping. We look into randomized sampling techniques which generate stationary samples. As an example, random walk based algorithms are in use in many cases because of several advantages offered by them [START_REF] Avrachenkov | Improving random walk estimation accuracy with uniform restarts[END_REF][START_REF] Brin | The anatomy of a large-scale hypertextual web search engine[END_REF].

We focus on the extremal properties in the correlated and stationary sequence of characteristics of interest X 1 , . . . , X n which is a function of the node sequence, the one actually generated by sampling algorithms. The characteristics of interest, for instance, can be node degrees, node income, number of followers of the node in OSN etc. Among the properties, clusters of exceedances of such sequences over high thresholds are studied in particular. The cluster of exceedances is determined as the consecutive exceedances of {X n } over the threshold {u n } between two consecutive non-exceedances [START_REF] Ferro | Inference for clusters of extreme values[END_REF][START_REF] Markovich | Modeling clusters of extreme values[END_REF]. It is important to investigate stochastic nature of extremes since it allows us to disseminate advertisement or collect opinions more effectively within the clusters.

The dependence structure of sampled sequence exceeding sufficiently high thresholds is measured using a parameter called extremal index (EI), θ in Extremal Value Theory. It is defined as follows.

Definition 1. [15, p. 53] The stationary sequence {X n } n≥1 , with F as the marginal distribution function and M n = max{X 1 , ..., X n }, is said to have the extremal index θ ∈ [0, 1] if for each 0 < τ < ∞ there is a sequence of real numbers (thresholds) u n = u n (τ ) such that

lim n→∞ n(1 -F (u n )) = τ and
(1)

lim n→∞ P{M n ≤ u n } = e -τ θ .
The maxima M n is related to EI more clearly as [6, p. 381] (a)

P{M n ≤ u n } = F nθ (u n ) + o(1). (2)
When {X n } n≥1 is i.i.d. (for instance uniform independent node sampling), θ = 1 and point processes of exceedances over threshold u n converges weakly to homogeneous Poisson process [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF]Chapter 5]. But when 0 ≤ θ < 1, point processes of exceedances converges weakly to compound Poisson process and this implies that exceedances of high threshold values u n tend to occur in clusters for dependent data [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF]Chapter 10]. EI has many useful interpretations and applications like

• Finding distribution of order statistics of the sampled sequence. These can be used to find quantiles and predicts the kth largest value which arise with a certain probability. Specifically for the distribution of maxima, (2) is available and the quantile of maxima is proportional to EI. Hence in case of samples with lower EI, lower values of maxima can be expected. When sampled sequence is the sequence of node degrees, these give many useful results.

• Close relation of extremal index to the distribution and expectation of the size of clusters of exceedances.

• First hitting time of the sampled sequence to (u n , ∞) is related to EI. Thus in case of applications where the aim is to detect large values of samples quickly, without actually employing sampling (which might be very costly), we can compare different sampling procedures by EI: smaller EI leads to longer searching of the first hitting time.

These interpretations are explained later in the paper. The network topology determines the stationary distribution of the characteristics of interest under a sampling technique and is reflected on the EI. This indicates that different sampling algorithms may have different EI.

Our contributions

The main contributions in this work are as follows. We associated Extremal Value Theory of stationary sequences to sampling of large complex networks and we study the extremal and clustering properties of the sampling process due to correlations. In order to facilitate a painless future study of correlations and clusters of samples in large networks, we propose to abstract the extremal properties into a single and handy parameter, EI. For any general stationary samples meeting two mixing conditions, we find that knowledge of bivariate distribution or bivariate copula is sufficient to compute EI analytically and thereby deriving many extremal properties. Several useful applications of EI (first hitting time, order statistics and mean cluster size) to analyse large graphs, known only through sampled sequences, are proposed. Degree correlations are explained in detail with a random graph model for which joint degree correlations exist for neighbor nodes. Three different random walk based algorithms that are widely discussed in literature (see [START_REF] Avrachenkov | Improving random walk estimation accuracy with uniform restarts[END_REF] and the references therein), are then revised for degree state space and EI is calculated when the joint degree correlation is bivariate Pareto distributed. We establish a general lower bound for EI in PageRank processes irrespective of the degree correlation model. Finally two estimation techniques of EI are provided and EI is numerically computed for a synthetic graph with neighbour degrees correlated and for two real networks (Enron email network and DBLP network).

The paper is organized as follows. In Section 2, methods to derive EI are presented. Section 3 considers the case of degree correlations. In Section 3.1 the graph model and correlated graph generation technique are presented. Section 3.2 explains the different types of random walks studied and derives associated transition kernels and joint degree distributions. EI is calculated for different sampling techniques later in Section 3.3. In Section 4 we provide several applications of extremal index in graph sampling techniques. In Section 5 we estimate extremal index and perform numerical comparisons. Finally Section 6 concludes the paper.

A shorter version of this submission has been appeared in [START_REF] Avrachenkov | Distribution and dependence of extremes in network sampling processes[END_REF].

Calculation of Extremal Index (EI)

We consider networks represented by an undirected graph G with N vertices and M edges. Since the networks under consideration are huge, we assume it is impossible to describe them completely, i.e., no adjacency matrix beforehand. Assume any randomized sampling procedure is employed and let the sampled sequence {X i } be any general sequence. This section explains a way to calculate extremal index from the bivariate distribution if the sampled sequence admits two mixing conditions.

Inria

Condition (D(u n )).

P(X i1 ≤ u n , . . . , X ip ≤ u n , X j1 ≤ u n , . . . , X jq ≤ u n ) -P(X i1 ≤ u n , . . . , X ip ≤ u n )P(X j1 ≤ u n , . . . , X jq ≤ u n ) ≤ α n,ln ,
where α n,ln → 0 for some sequence l n = o(n) as n → ∞, for any integers i 1 ≤ . . . < i p < j 1 < . . . ≤ j q with j 1 -i p > l n .

Condition (D ′′ (u n )). lim n→∞ rn j=2 P(X j ≤ u n < X j+1 |X 1 > u n ) = 0,
where (n/r n )α n,ln → 0 and l n /r n → 0 with α n,ln , l n as in Condition D(u n ) and r n as o(n).

Let C(u, v) is the bivariate Copula [19] ([0, 1] 2 → [0, 1]
) and C ′ is its Gâteaux derivative along the direction (1, 1). Using Sklar's theorem [19, p. 18], with F as the marginal stationary distribution function of the sampling process,

C(u, u) = P(X 1 ≤ F -1 (u), X 2 ≤ F -1 (u)).
F -1 denotes the inverse function of F . This representation is unique if the stationary distribution F (x) is continuous.

Proposition 1. If the sampled sequence is stationary and satisfies conditions D(u n ) and D ′′ (u n ), then extremal index is given by θ

= C ′ (1, 1) -1, (3)
and 0 ≤ θ ≤ 1.

Proof. From [START_REF] Leadbetter | On exceedance point processes for stationary sequences under mild oscillation restrictions[END_REF], for the stationary sequence {X n } with Conditions D(u n and

D ′′ (u n ), θ = lim n→∞ P(X 2 ≤ u n |X 1 > u n ). Then θ = lim n→∞ P(X 2 ≤ u n , X 1 > u n ) P(X 1 > u n ) = lim n→∞ P(X 2 ≤ u n ) -P(X 1 ≤ u n , X 2 ≤ u n ) P(X 1 > u n ) = lim n→∞ P(X 2 ≤ u n ) -C P(X 1 ≤ u n ), P(X 2 ≤ u n ) 1 -P(X 1 ≤ u n ) = lim x→1 x -C(x, x) 1 -x = C ′ (1, 1) -1.
The existence of EI in [0, 1] is evident from the definition used in this proof.

Remark 1. The condition D ′′ (u n ) can be made weaker to D (k) (u n ) presented in [START_REF] Chernick | Calculating the extremal index for a class of stationary sequences[END_REF],

lim n→∞ nP X 1 > u n ≥ max 2≤i≤k X i , max k+1≤j≤rn X j > u n = 0, RR n°8578
where r n is defined as in D ′′ (u n ). For the stationary sequence D (2) k) is satisfied for some k ≥ 2 along with D(u n ), then following the proof of Proposition 1, EI can be derived as

(u n ) = D ′′ (u n ). If we assume D (
θ = C ′ k (1) -C ′ k-1 (1), where C k (x) represents the copula of k-dimensional vector (x 1 , . . . , x k ), C k (x 1 , . . . , x k ) with x = x 1 . . . = x k and C k-1 is its (k -1)th marginal, C k-1 (x) = C k-1 (x 1 , . . . , x k-1 , 1) with x = x 1 . . . = x k-1 .
In some cases it is easy to handle with the joint tail distribution. Survival Copula C(•, •) which corresponds to

P(X 1 > x, X 2 > x) = C(F (x), F (x)), with F (x) = 1 -F (x), can also be used to calculate θ. It is related to Copula as C(u, u) = C(1 -u, 1 -u) + 2u -1 [19, p. 32]. Hence θ = C ′ (1, 1) -1 = 1 -C ′ (0, 0).
Lower tail dependence function of survival copula is defined as [START_REF] Weng | Characterization of multivariate heavy-tailed distribution families via copula[END_REF] λ(u 1 , u 2 ) = lim

t→0 + C(tu 1 , tu 2 ) t .
Hence C ′ (0, 0) = λ(1, 1). λ can be calculated for different copula families. In particular, if C is a bivariate Archimedean copula, then it can be represented as,

C(u 1 , u 2 ) = ψ(ψ -1 (u 1 ) + ψ -1 (u 2 )),
where ψ is the generator function and

ψ -1 is its inverse with ψ : [0, ∞] → [0, 1] meeting several other conditions. If ψ is a regularly varying distribution with index -β, β > 0, then λ(x 1 , x 2 ) = (x -β -1 1 +x -β -1 2
) -β and (X 1 , X 2 ) has a multivariate regularly varying distribution [START_REF] Weng | Characterization of multivariate heavy-tailed distribution families via copula[END_REF]. Therefore, for Archimedean copula family, EI is given by

θ = 1 -1/2 β . (4)
As an example, bivariate Pareto distribution of the form P(

X 1 > x 1 , X 2 > x 2 ) = (1 + x 1 + x 2 ) -γ , γ > 0 has Arhimedean copula with generator function ψ(x) = (1 + x) -γ . This gives θ = 1 -1/2 γ .
Bivariate exponential distribution of the form

P(X 1 > x 1 , X 2 > x 2 ) = 1 -e -x1 -e -x2 + e -(x1+x2+ηx1x2) , 0 ≤ η ≤ 1, also admits Archimedian copula.
2.1 Check of conditions D(u n ) and D ′′ (u n )

If the sampling technique is assumed to be based on a Markov chain and consider the sampled sequence as measurable functions of stationary Markov samples, then such a sequence is stationary and [START_REF] O'brien | Extreme values for stationary and markov sequences[END_REF] proved that another mixing condition AIM (u n ) which implies D(u n ) is satisfied. Condition D ′′ (u n ) allows clusters with consecutive exceedances and eliminates the possibility of clusters with upcrossing of the threshold u n (X i ≤ u n < X i+1 ). Hence in those cases, where it is tedious to check the condition D ′′ (u n ) theoretically, we can use numerical procedures to measure ratio of number of consecutive exceedances to number of exceedances and the ratio of number of upcrossings to number of consecutive exceedances in small intervals. Such an example is provided in Section 3.3.

Remark 2. The EI is derived in [START_REF] Ferreira | Extremal functions, extremal index and markov chains[END_REF] to the same expression in [START_REF] Avrachenkov | Distribution and dependence of extremes in network sampling processes[END_REF]. But [START_REF] Ferreira | Extremal functions, extremal index and markov chains[END_REF] assumes {X n } is sampled from a first order Markov chain. This condition is much stricter than D(u n ) and D ′′ (u n ) which we used to derive (3). For instance, degrees of the node samples obtained from Inria a Markov chain based sampling, mostly not form a Markov chain as node-degree relation is not one-one while D(u n ) is agreed for such a case and D ′′ (u n ) can get satisfied, see Section 3.3 for an example.

Degree correlations

The techniques established in Section 2 are very general, applicable to any sampling techniques and any sequence of samples which satisfy certain conditions. In this section we illustrate the calculation of extremal index for correlations among degrees. We introduce different sampling techniques through this section though they can be used in case of any general correlations. We denote the sampled sequence {X i } as {D i } in this section.

Description of the model

We take into account correlation in degrees between neighbor nodes. The dependence structure in the graph is described by the joint degree-degree probability density function f (d 1 , d 2 ) with d 1 and d 2 indicating the degrees of adjacent nodes or equivalently by the corresponding tail distribution function

F (d 1 , d 2 ) = P(D 1 ≥ d 1 , D 2 ≥ d 2 )
with D 1 and D 2 representing the corresponding degree random variables (see e.g., [START_REF] Barrat | Dynamical processes on Complex Networks[END_REF][START_REF] Boguna | Epidemic spreading in complex networks with degree correlations[END_REF][START_REF] Goltsev | Percolation on correlated networks[END_REF]).

The probability that a randomly chosen edge has the end vertices with degrees

d 1 ≤ d ≤ d 1 + ∆(d 1 ) and d 2 ≤ d ≤ d 2 + ∆(d 2 ) is (2 -δ d1d2 )f (d 1 , d 2 )∆(d 1 )∆(d 2 ). Here δ d1d2 = 1 if d 1 = d 2 ,
zero otherwise. The multiplying factor 2 appear on the above expression when

d 1 = d 2 because of the symmetry in f (d 1 , d 2 ), f (d 1 , d 2 ) = f (d 2 , d 1 )
due to the undirected nature of the underlying graph, and the fact that both f (d 1 , d 2 and f (d 2 , d 1 ) contribute to the edge probability under consideration.

The degree density f d (d 1 ) can be calculated from the marginal of f (d 1 , d 2 ) as

f (d 1 ) = d2 f (d 1 , d 2 )d(d 2 ) ≈ d 1 f d (d 1 ) E[D] , (5)
where E[D] denotes the mean node degree,

E[D] = f (d 1 , d 2 ) d 1 d(d 1 )d(d 2 ) -1
. f (.) can be interpreted as the degree density of a vertex reached by following a randomly chosen edge. The approximation for f (d 1 ) is obtained as follows: in the R.H.S. of (5), roughly,

d 1 f d (d 1 )
N is the number of half edges from nodes with degree around d 1 and E[D]N is the total number of half edges.

From the above description, it can be noted that the knowledge of f (d 1 , d 2 ) is sufficient to describe this random graph model and for its generation.

Most of the results in this paper are derived assuming continuous probability distributions for f (d 1 , d 2 ) and f d (d 1 ) because an easy and unique way to calculate extremal index exists for continuous distributions in our setup (more details in Section 2). Also the extremal index might not exist for many discrete valued distributions [START_REF] Leadbetter | Extremes and related properties of random sequences and processes[END_REF].

Random graph generation

A random graph bivariate joint degree-degree correlation distribution can be generated as follows ([20]).

1. Degree sequence is generated according to the degree distribution,

f d (d) = f (d)E[D] d 2.
An uncorrelated random graph is generated with the generated degree sequence using configuration model ([5])

3. Metropolis dynamics is applied now on the generated graph: choose two edges randomly (denoted by the vertex pairs (v 1 , w 1 ) and (v 2 , w 2 )) and measure the degrees, (j 1 , k 1 ) and (j 2 , k 2 ) correspond to these vertex pairs. Generated a random number, y, according to uniform distribution in

[0, 1]. If y ≤ min(1, (f (j 1 , j 2 )f (k 1 , k 2 ))/(f (j 1 , k 1 )f (j 2 , k 2 )
)), then remove the selected edges and construct news ones as (v 1 , v 2 ) and (w 1 , w 2 ). Otherwise keep the selected edges intact. This dynamics will generate the required joint degree-degree distribution. Run Metropolis dynamics well enough to mix the network.

Description of random walks

In this section, we explain three different random walk based algorithms for exploring the network. They have been extensively studied in previous works [START_REF] Avrachenkov | Improving random walk estimation accuracy with uniform restarts[END_REF][START_REF] Brin | The anatomy of a large-scale hypertextual web search engine[END_REF][START_REF] Lovász | Random walks on graphs: A survey[END_REF] where they are formulated with vertex set as the state space of the underlying Markov chain on graph. The walker in these algorithms, after reaching each node, moves to another node randomly by following the transition kernel of the Markov chain. But since the interest in the present work is in the degree sequence, rather than node sequence, and its extremal properties, we take degree set as the state space and find appropriate transition kernels. We use f X and P X to represent the probability density function and probability measure under the algorithm X with the exception that f d represents the probability density function of degrees.

Random Walk (RW)

In a random walk, the next node to visit is chosen uniformly among the neighbors of the current node. From (5) we approximate the standard random walk on degree state space by the following transition kernel, conditional density function that the present node has degree d t and the next node is with degree d t+1 ,

f RW (d t+1 |d t ) ≈ E[D]f (d t , d t+1 ) d t f d (d t ) . ( 6
)
This approximation is obtained as follows: given the present node has degree d t , 1/d t is the probability of selecting a neighbor uniformly and rest of the terms in R.H.S. represent the mean number of neighbors with degree around d t+1 . When

d t = d t+1 , E[D]N 2 (2f (d t , d t+1 ))
is the mean number of edges between degrees about d t and d t+1 and f d (d t )N is the mean number of nodes with degrees about d t , and thus their ratio represents such a mean number of edges per node with degree about d t , i.e., mean number of neighbors with degree about d t+1 . The probability of occurring the other case, d t = d t+1 , is zero as the degrees are assumed to follow a continuous distribution.

If the standard random walk on the vertex set is in the stationary regime, its stationary distribution (probability of staying at a particular vertex i) is proportional to the degree (see e.g., [START_REF] Lovász | Random walks on graphs: A survey[END_REF]) and is given by d i /2M . Then in the standard random walk on degree set, the stationary distribution of staying at any node with degree around d 1 can be approximated as

N f d (d 1 ) (d 1 /2M |). Thus f RW (d 1 ) ≈ d 1 E[D] f d (d 1 ),
Then, the joint density of the standard random walk is

f RW (d t+1 , d t ) ≈ f (d t , d t+1 ).

Inria

Check of the approximation

We provide comparison of simulated values and theoretical values of transition kernel of RW in Figure 1. The bivariate Pareto model is assumed for the joint degree-degree tail function of the graph, PageRank is a modification of the random walk which with a fixed probability 1 -c samples a random node with uniform distribution and with a probability c, it follows the random walk transition [START_REF] Brin | The anatomy of a large-scale hypertextual web search engine[END_REF]. Its evolution on degree state space can be described as follows:

F (d 1 , d 2 ) = 1 + d 1 -µ σ + d 2 -µ σ -γ , (7)
f P R (d t+1 |d t ) ≈ c f RW (d t+1 |d t ) + (1 -c) 1 N N f d (d t+1 ) ≈ c f RW (d t+1 |d t ) + (1 -c)f d (d t+1 ) (8)
Here the 1/N corresponds to the uniform sampling on vertex set and 1 N N f d (d t+1 ) indicates the net probability of jumping to all the nodes with degree around d t+1 .

Check of the approximation

We provide a consistency check of the approximation derived for transition kernel by studying tail behavior of degree distribution and PageRank distribution. It is known that under some strict conditions, for a directed graph, PageRank and Indegree have same tail exponents [START_REF] Litvak | In-degree and pagerank: Why do they follow similar power laws?[END_REF]. In our formulation in terms of degrees, for uncorrelated and undirected graph, PageRank for a given degree d, P R(d), can be approximated from the basic definition as,

P R(d) = f RW (d) ≈ c f RW (d) + (1 -c) f d (d).

RR n°8578

This is a deterministic quantity. We are interested in the distribution of the random variable P R(D), PageRank of a randomly choosen degree class D. PageRank P R(d) is also the long term proportion or probability that PageRank process ends in a degree class with degree d. This can be scaled suitably to provide a rank-type information. Its tail distribution is

P (P R(D) > x) = P (c.f RW (D) + (1 -c).f d (D) > x) ,
where D ∼ f d (.). The PageRank of any vertex inside the degree class d is P R(d)/(N f d (d)). The distribution of Page Rank of a randomly chosen vertex i, P (P R(i) > x) after appropriate scaling for comparison with degree distribution is P (N.P R(i) > d), where d = N x. Now

P (N.P R(i) > d) = P N P R(D) N f d (D) > d = P D > E[D] c d -(1 -c) .
This of the form P (D > A d + B) with A and B as appropriate constants and hence will have the same exponent of degree distribution tail when the graph is uncorrelated.

There is no convenient expression for the stationary distribution of PageRank, to the best of our knowledge, and it is difficult to come up with an easy to handle expression for the joint distribution. Therefore, along with other advantages, we consider another modification of the standard random walk.

Random Walk with Jumps (RWJ)

RW sampling leads to many practical issues like the possibility to get stuck in a disconnected component, biased estimators etc. RWJ overcomes such problems ( [START_REF] Avrachenkov | Improving random walk estimation accuracy with uniform restarts[END_REF]).

In this algorithm we follow random walk on a modified graph which is a superposition of the given graph and complete graph on same vertex set of the given graph with weight α/N on each edge, α ∈ [0, ∞] being a design parameter ( [START_REF] Avrachenkov | Improving random walk estimation accuracy with uniform restarts[END_REF]). The algorithm can be shown to be equivalent to select c = α/(d t + α) in the PageRank algorithm, where d t is the degree of the present node. The larger the node's degree, less likely is the artificial jump of the process. This modification makes the underlying Markov chain time reversible, significantly reduces mixing time, improves estimation error and leads to a closed form expression for stationary distribution.

The transition kernel on degree set, following PageRank kernel, is

f RW J (d t+1 |d t ) ≈ d t d t + α f RW (d t+1 |d t ) + α d t + α f d (d t+1 ) = E[D]f (d t , d t+1 ) + αf d (d t )f d (d t+1 ) (d t + α)f d (d t ) .
The stationary distribution for node i (on the vertex set) is (d i + α)/(2M + N α) and the equivalent stationary probability density function on degree set by collecting all the nodes with same degree is

f RW J (d 1 ) ≈ d 1 + α 2M + N α N f d (d 1 ) = (d 1 + α)f d (d 1 ) E[D] + α , (9)
Inria since 2M/N = E[D]
. The stationarity of the f RW J (d 1 can be verified by plugging the obtained expression in the stationarity condition of the Markov Chains. We have

f RW J (d 1 ) = f RW J (d 1 |d 2 )f RW J (d 2 )d(d 2 ) ≈ E[D]f (d 1 , d 2 ) + αf d (d 1 )f d (d 2 ) (d 2 + α)f d (d 2 ) (d 2 + α)f d (d 2 ) E[D] + α d(d 2 ) ≈ (d 1 + α)f d (d 1 ) E[D] + α ,
where (5) has been applied. Then, the joint density function for the random walk with jumps has the following form

f RW J (d t+1 , d t ) ≈ E[D]f (d t+1 , d t ) + αf d (d t+1 )f d (d t ) E[D] + α .
Moreover the associated tail distribution has a simple form,

f RW J (D t+1 > d t+1 , D t > d t ) ≈ E[D]F (d t+1 , d t ) + αF d (d t+1 )F d (d t ) E[D] + α . ( 10
)
Remark 3. Characterizing Markov chain based sampling in terms of degree transition has some advantages,

• In the different random walk algorithms considered on vertex set, all the nodes with same degree have same stationary distribution. This also implies that it is more natural to formulate the random walk transition in terms of degree.

• Degree uncorrelations in the underlying graph is directly reflected in the joint distribution of the studied sampling techniques. For uncorrelated networks,

f RW (d 1 , d 2 ) = f RW (d 1 )f RW (d 2 ), f P R (d 1 , d 2 ) = f P R (d 1 )f P R (d 2 ) and f RW J (d 1 , d 2 ) = f RW J (d 1 )f RW J (d 2 ).

Extremal Index for bivariate Pareto Degree Correlation

As explained in the Introduction section, extremal index is an important parameter in characterizing dependence and extremal properties in a stationary sequence. We assume that we have waited sufficiently long that the underlying Markov chain of the three different graph sampling algorithms are in stationary regime now. Here we derive EI of RW and RWJ for the model with degree correlation among neighbours as bivariate Pareto [START_REF] Boguna | Epidemic spreading in complex networks with degree correlations[END_REF]. The two mixing conditions D(u n ) and D ′′ (u n ) introduced in Section 2 are needed for our EI analysis. Condition D(u n ) is satisfied as explained in Section 2.1. An empirical evaluation of D ′′ (u n ) is provided in Section 5.3.1.

EI for Random Walk sampling

We use the expression for EI given in Proposition 1. As f RW (x, y) is same as f (x, y), we have, γ+1) .

C(u, u) = P(D 1 > F -1 (u), D 2 > F -1 (u)) = 1 + 2(u -1/γ -1) -γ C ′ (u, u) = 2(2 -u 1/γ ) -(
Thus θ = 1 -C ′ (0, 0) = 1 -1/2 γ . For γ = 1 we get θ = 1/2.
In this case, we can also use expression given in (4).

EI for Random Walk with Jumps sampling

Although it is possible to derive EI as in RW case above, we provide an alternative way to avoid the calculation of tail distribution of degrees and inverse of RWJ marginal (with respect to the bivariate Pareto degree correlation).

Under the assumption of D ′′ ,

θ = lim n→∞ P (D 2 ≤ u n , D 1 > u n ) P (D 1 > u n ) = lim n→∞ P (D 1 ≥ u n ) -P (D 2 ≥ u n , D 1 ≥ u n ) P (D 1 > u n ) (11)
Now using the condition (1) on the marginal and joint tail distribution of RWJ [START_REF] Ferreira | Extremal functions, extremal index and markov chains[END_REF], we can write (b)

P (D 1 ≥ u n ) -P (D 2 ≥ u n , D 1 ≥ u n ) P (D 1 > u n ) = τ /n + o(1/n) -E[D] E[D]+α P RW (D 2 ≥ u n , D 1 ≥ u n ) - α E[D]+α O(τ /n)O(τ /n) τ /n + o(1/n)
The asymptotics in the last term of the numerator is due to the following:

F RW J (u n ) = E[D] E[D] + α F (u n ) + α E[D] + α F d (u n ) = τ /n + o(1/n),
and hence

F d (u n ) = O(τ /n). Therefore (11) becomes θ = 1 - E[D] E[D] + α lim n→∞ P RW (D 2 ≥ u n , D 1 ≥ u n )n/τ
In the case of the bivariate Pareto distribution [START_REF] Boguna | Epidemic spreading in complex networks with degree correlations[END_REF], we obtain

θ = 1 - E[D] E[D] + α 2 -γ (12)

Lower bound of EI of the PageRank

We obtain the following lower bound for EI in the PageRank processes.

Proposition 2. For the PageRank process on degree state space irrespective of the degree correlation structure in the underlying graph, the extremal index

θ ≥ (1 -c).
Proof. From [START_REF] O'brien | Extreme values for stationary and markov sequences[END_REF], the following representation of EI holds for degree sequence,

lim n→∞ P{M 1,pn ≤ u n |D 1 > u n } = θ, ( 13
)
where {p n } is an increasing sequence of positive integers, p n = o(n) as n → ∞ and M 1,pn = max{D 2 , ..., D pn }. Let A be the event that the node corresponding to D 2 is selected uniformly Inria among all the nodes, not following random walk from the node for D 1 . Then P P R (A) = 1 -c. Now, with (8), 1)

P P R (M 1,pn ≤ u n |D 1 > u n ) ≥ P P R (M 1,pn ≤ u n , A|D 1 > u n ) = P P R (A|D 1 > u n )P P R (M 1,pn ≤ u n |A, D 1 > u n ) (i) = (1 -c)P P R (M 1,pn ≤ u n ), (ii) = (1 -c)P (pn-1)θ P R (D 1 ≤ u n ) + o(
≥ (1 -c)P (pn-1) P R (D 1 ≤ u n ) + o(1) (iii) ∼ (1 -c)(1 -τ /n) pn-1 , ( 14
)
where {p n } is the same sequence as in ( 13) and (i) follows mainly from the observation that conditioned on A, {M 1,pn ≤ u n } is independent of {D 1 > u n }, (ii) and (iii) result from the approximations in ( 2) and ( 1) respectively. Assuming p n -1 = n 1/2 and since (1 13) and ( 14), θ ≥ 1 -c.

-τ /n) pn-1 ∼ e -τ / √ n → 1 as n → ∞, from (
The PageRank transition kernel (8) on the degree state space does not depend upon the random graph model in Section 3.1. Hence the derived lower bound of EI is useful for any degree correlation model.

Applications of Extremal Index in Network Sampling Processes

This section provides several uses of EI to infer the sampled sequence. This emphasis that the analytical calculation and estimation of EI are practically relevant. The limit of the point process of exceedances, N n (.), which counts the times, normalized by n, at which {X i } n i=1 exceeds a threshold u n provides many applications of extremal index. A cluster is considered to be formed by the exceedances in a block of size r n (r n = o(n)) in n with cluster size ξ n = rn i=1 1(X i > u n ) when there is at least one exceedance within r n . The point process N n converges weakly to a compound poisson process (CP ) with rate θτ and i.i.d. distribution as the limiting distribution of cluster size, under condition (1) and a mixing condition, and the points of exceedances in CP correspond to the clusters [6, Section 10.3]. We name this kind of clusters as blocks of exceedances.

The applications below require a choice of the threshold sequence {u n } satisfying (1). For practical purposes, if a single threshold u is demanded for the sampling budget B, we can fix u = max{u 1 , . . . , u B }.

The applications in this section are explained with the assumption that the sampled sequence is the sequence of node degrees. But the following techniques are very general and can be extended to any sampled sequence satisfying conditions D(u n ) and D ′′ (u n ).

Order statistics of the sampled degrees

The order statistics X n-k,n , (n -k)th maxima, is related to N n (.) and thus to θ by

P(X n-k,n ≤ u n ) = P(N n ((0, 1]) ≤ k),
where we apply the result of convergence of N n to CP [6, Section 10.3.1].

Distribution of Maxima

The distribution of the maxima of the sampled degree sequences can be derived as (2) when n → ∞. Hence if the extremal index of the underlying process is known then from (2) one can approximate the (1 -η)th quantile x η of the maximal degree M n as

P{M n ≤ x η } = F nθ (x η ) = P nθ {X 1 ≤ x η } = 1 -η, i.e. x η ≈ F -1 (1 -η) 1/(nθ) . (15)
In other words, quantiles can be used to find the maxima of the degree sequence with certain probability.

For a fixed certainty η, x η is proportional to θ. Hence if the sampling procedures have same marginal distribution, with calculation of EI, it is possible to predict how much large values can be achieved. Lower EI indicates lower value for x η and higher represents high x η . For the random walk example in Section 3.3.1 for the degree correlation model, with the use of ( 15), we get the (1 -η)th quantile of the maxima M n

x η ≈ µ + σ 1 -(1 -η) 1/(nθ) -1/γ -1
The following example demonstrates the effect of neglecting correlations on the prediction of the largest degree node. The largest degree, with the assumption of Pareto distribution for the degree distribution, can be approximated as KN 1/δ with K ≈ 1, N as the number of nodes and γ as the tail index of complementary distribution function of degrees [START_REF] Avrachenkov | Quick detection of nodes with large degrees[END_REF]. For Twitter graph (recorded in 2012), δ = 1.124 for outdegree distribution and N = 537, 523, 432 [START_REF] Gabielkov | Studying social networks at scale: Macroscopic anatomy of the twitter social graph[END_REF]. This gives the largest degree prediction as 59, 453, 030. But the actual largest out degree is 22, 717, 037. This difference is because the analysis in [START_REF] Avrachenkov | Quick detection of nodes with large degrees[END_REF] assumes i.i.d. samples and does not take into account the degree correlation. With the knowledge of EI, correlation can considered as in (2). In the following section, we derive an expression for such a case.

Estimation of largest degree when the marginals are Pareto distributed

It is known that many social networks have the degree asymptotically distributed as Pareto. We find that in these cases, the marginal distribution of degrees of the random walk based methods also follow Pareto distribution (though we have derived only for the model with degree correlations among neighbors, see Section 3) Claim 1. For any stationary sequence with marginal distribution following Pareto distribution F (x) = Cx -δ , the largest value is

M n ≈ (nθ) 1/δ C log 2 1/δ
Proof. From extreme value theory [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF], it is known that when {X i , i ≥ 1} are i.i.d.,

lim n→∞ P M n -b n a n ≤ x = H γ (x), (16)
where H γ (x) is the extreme value distribution with index γ and {a n } and {b n } are appropriately chosen deterministic sequences. When {X i , i ≥ 1} are stationary with EI θ, the limiting distribution becomes H ′ γ ′ (x) and it differs from H γ (x) only through parameters. H γ (x) = exp(-t(x))

Inria with t(x) = 1 + x-µ σ γ -1/γ
. With the normalizing constants (µ = 0 and σ = 1), H ′ γ ′ has the same shape as H γ with parameters γ ′ = γ, σ ′ = θ γ and µ ′ = (θ γ -1)/γ.

For Pareto case, F (x) = Cx -δ , γ = 1/δ, a n = γC γ n γ and b n = C γ n γ . From ( 16), for large n, M n is stochastically equivalent to a n χ + b n , where χ is a random variable with distribution H ′ γ ′ . It is observed in [START_REF] Avrachenkov | Quick detection of nodes with large degrees[END_REF] that median of χ is an appropriate choice for the estimation of

M n . Median of χ = µ ′ + σ ′ (log 2) -γ ′ -1 γ ′ = (θ γ (log 2) -γ -1)γ -1 . Hence, M n ≈ a n θ γ (log 2) -γ γ -1 + b n = (nθ) 1/δ C log 2 1/δ

Relation to first hitting time and interpretations

Extremal index also gives information about the first time {X n } hits (u n , ∞). Let T n be this time epoch. As N n converges to compound poisson process, it can be observed that T n /n is asymptotically an exponential random variable with rate θτ , i.e., lim n→∞ P(T n /n > x) = exp(-θτ x). Therefore lim n→∞ E(T n /n) = 1/(θτ ). Thus the more EI smaller, the more time it will take to hit the extreme levels as compared to independent sampling. This property can make use to compare different sampling procedures.

Relation to mean cluster size

If the conditions D ′′ (u n ) is satisfied along with D(u n ), asymptotically, a run of the consecutive exceedances following an upcrossing is observed, i.e., {X n } crosses the threshold u n at a time epoch and stays above u n for some more time before crossing u n downwards and stays below it for some time until next upcrossing of u n happens. This is called cluster of exceedances and is more practically relevant than blocks of exceedances at the starting of this section and is shown in [START_REF] Leadbetter | On exceedance point processes for stationary sequences under mild oscillation restrictions[END_REF] that these two definitions clusters are asymptotically equivalent resulting in similar cluster size distribution. The expected value of cluster of exceedances converges to inverse of extremal index [6, p. 384], i.e.,

θ -1 = lim n→∞ j≥1
jπ n (j),

where {π n (j), j ≥ 1} is the distribution of size of cluster of exceedances with n samples. More details about cluster size distribution and its mean can be found in [START_REF] Markovich | Modeling clusters of extreme values[END_REF].

Estimation of Extremal Index and Numerical results

This section introduces two estimators for EI. Two types of networks are presented: synthetic correlated graph and real networks (Enron email network and DBLP network). For the synthetic graph, we compare the estimated EI to its theoretical value. For the real network, we calculate EI using the two estimators. We take {X i } as the degree sequence and use RW, PR and RWJ as the sampling techniques. The methods mentioned in the following are general and are not specific to degree sequence or random walk technique.

Empirical Copula based estimator

We have tried different estimators for EI available in literature [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF][START_REF] Ferreira | Extremal functions, extremal index and markov chains[END_REF] and found that the idea of estimating copula and then finding value of its derivative at (1, 1) works without the need to choose and optimize several parameters found in other estimators. We assume that {X i } satisfies D(u n ) and D ′′ (u n ) and we use (3) for calculation of EI. Copula C(u, v) is estimated empirically by

C n (u, v) = 1 n n k=1 I R X i k n + 1 ≤ u, R Y i k n + 1 ≤ v , with R X i k indicates rank of the element X i k in {X i k , 1 ≤ k ≤ n}, and Y i k = X i k +1 .
The sequence {X i k } is chosen from the original sequence {X i } in such a way that X i k and X i k+1 are sufficiently apart to make them independent to certain extent. The large-sample distribution of C n (u, v) is normal and centered at copula C(u, v). Now, to get θ, we use linear least squares error fitting to find slope at (1, 1) or use cubic spline interpolation for better results.

Intervals Estimator

This estimator does not assume any conditions on {X i }, but has the parameter u to choose appropriately. Let N = n i=1 1(X i > u) be number of exceedances of u at time epochs 1 ≤ S 1 < . . . < S N ≤ n and let the interexceedance times are T i = S i+1 -S i . Then intervals estimator is defined as [6, p. 391],

θn (u) = min(1, θ1 n (u)), if max T i : 1 ≤ i ≤ N -1 ≤ 2, min(1, θ2 n (u)), if max T i : 1 ≤ i ≤ N -1 > 2,
where

θ1 n (u) = 2( N -1 i=1 T i ) 2 (N -1) N -1 i=1 T 2 i , and
θ2 n (u) = 2( N -1 i=1 (T i -1)) 2 (N -1) N -1 i=1 (T i -1)(T i -2)
.

We choose u as δ percentage quantile thresholds, i.e., δ percentage of {X i , 1 ≤ i ≤ n} falls below u. The EI is usually selected corresponding to the stability interval in the plot (θ, δ).

Synthetic graph

The simulations in the section follow the bivariate Pareto model and parameters introduced in (7). We use the same set of parameters of Figure 1 and the graph is generated according to the technique in Section 3.1.1.

For the RW case, Figure 2 shows copula estimator, and theoretical copula based on the continuous distribution in [START_REF] Boguna | Epidemic spreading in complex networks with degree correlations[END_REF] and is given by

C(u, u) = 1 + 2((1 -u) -1/γ -1) -γ + 2u -1.
Though we take quantized values for degree sequence, it is found that the copula estimated matches with theoretical copula. The value of EI is then obtained after cubic interpolation and numerical differentiation at point (1, 1). For the theoretical copula, EI is 1 -1/2 γ , where γ = 1.2. For the RWJ algorithm, Figure 4 shows the Intervals estimate and theoretical value for different α. We used the expression [START_REF] Gabielkov | Studying social networks at scale: Macroscopic anatomy of the twitter social graph[END_REF] for theoretical calculation. The small difference in theory and simulation results is due to the assumption of continuous degrees in the analysis, but the practical usage requires quantized version.

Figure 5 displays the Intervals estimate of EI with PR sampling. It can be seen that the lower bound proposed in Proposition 2 gets tighter as c decreases.

Check of condition D ′′

The mixing conditions D(u n ) and D ′′ (u n ) need to be satisfied for using the theory in Section 2. Though Intervals estimator does not require them, these conditions will provide the existence of EI. Condition D(u n ) works in this case as explained in previous sections and for D ′′ (u n ), we do the following empirical test. We collect samples for each of the techniques RW, PR and RWJ. Intervals are taken of duration 5, 10, 15 and 20 time samples. The ratio of number of upcrossings to number of exceedances r up and ratio of number consecuitve exceedances to number of exceedances r cluster are calculated in Table 1. These proportions are averaged over 2000 occurrences of each of these intervals and over all the different intervals. The statistics in the table indicates strong occurrence of condition D ′′ (u n ). We have also observed that the changes in the parameters does not affect this inference. Table 1: Test of Condition D ′′ in the synthetic graph

Real network

We consider two real world networks: Enron email network and DBLP network. The data is collected from [START_REF]Stanford large network dataset collection[END_REF]. Both the networks satisfy the check for the condition D ′′ (u n ) reasonably well.

For the RW sampling, Figure 6 shows the bivariate copula estimated and mentions corresponding EI. Intervals estimator is presented in Figure 7. After observing plateaus in the plots, we took EI as 0.25 and 0.2 for DBLP and Enron email graphs, respectively. lates this relation. We relate EI to many relevant extremes in networks like order statistics, first hitting time, mean cluster size etc. In particular, we model correlation in degrees of adjacent nodes and examine samples from random walks on degree state space. Finally we have obtained estimates of EI for a synthetic graph with degree correlations and find a good match with the theory. We also calculate EI for two real-world networks. In future, we plan to investigate the relation between assortativity coefficient and EI, and intends to study in detail the EI in real networks.

Endnotes

(a) F k (.) indicates kth power of F (.) throughout the paper except when k = -1 where it denotes the inverse function.

(b) ∼' stands for asymptotically equal, i.e. f (x) ∼ g(x) ⇔ f (x)/g(x) → 1 as x → a, x ∈ M where the functions f (x) and g(x) are defined on some set M and a is a limit point of M . f (x) = o(g(x)) means limx→a f (x)/g(x) = 0. Also f (x) = O(g(x)) indicates that there exist δ > 0 and M > 0 such that |f (x)| ≤ M |g(x)| for |x -a| < δ.

  where σ, µ and γ are positive values. In the figure, N number of nodes is 5,000. µ = 10, γ = 1.2 and σ = 15. These choices of parameters provides E[D] = 21.0052. At each instant Metropolis dynamics will choose two edges and it has run 200,000 times (provides sufficient mixing). The figure shows satisfactory fitting of the approximation.

Figure 1 :

 1 Figure 1: Transition kernel comparison

Figure 3

 3 displays the comparison between theoretical value of EI and Intervals estimate.

Figure 2 :

 2 Figure 2: Synthetic graph (RW sampling): Empirical and theoretical copulas

Figure 3 :

 3 Figure 3: Synthetic graph (RW sampling): Intervals estimate and theoretical value θ = 0.56 vs the percentage of quantile level, δ

Figure 4 :

 4 Figure 4: Synthetic graph (RWJ sampling): Intervals estimate and theoretical value vs the percentage of quantile level δ

Figure 5 :

 5 Figure 5: Synthetic graph (PR sampling): Intervals estimate vs the percentage of quantile level δ

Inria

  : θest. = 0.21 DBLP: θest. = 0.29

Figure 6 :

 6 Figure 6: Empirical copulas for email-Enron graph and DBLP graph (RW sampling)

Figure 7 :

 7 Figure 7: email-Enron graph and DBLP graph (RW sampling): Intervals estimate vs the percentage of quantile level δ.

  .

Figure 8 :

 8 Figure 8: email-Enron graph (RWJ sampling): Intervals estimate vs the percentage of quantile level δ.

  Distribution et Dépendance des Valeurs Extrêmes dans l'Échantillonnage des Réseaux Complexes Résumé : Nous étudions la structure de dépendance dans les séquences échantillonnés dans les réseaux complexes. Nous considérons des algorithmes randomisés pour échantillonner les noeuds et étudions les propriétés extrémales dans n'importe quelle suite des caractéristiques comme degrés de noeuds, nombres d'adeptes ou des revenus des noeuds des réseaux sociaux. Nous faisons l'abstraction de la dépendance des extrêmes par un seul paramètre qui apparaît dans la théorie des valeurs extrêmes. Ce paramètre s'appelle l'Indice Extrémal (IE). Dans ce travail, nous étudions ce paramètre analytiquement en utilisant copulas, et empiriquement par estimation. Nous proposons d'utiliser IE comme un paramètre pour comparer des différentes procédures d'échantillonnage. Comme un exemple, les corrélations de degré entre les noeuds voisins sont étudiées en détail en utilisant les trois techniques d'échantillonnage basé sur des marches aléatoires.

RR n°8578

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

The work of the first and third authors is partly supported by ADR "Network Science" of Alcatel-Lucent Inria joint lab. The second author is partly supported by the Russian Foundation for Basic Research, Grant 13-08-00744 A, and Campus France -Russian Embassy bilateral exchange program.

Conclusions

In this work, we have associated Extreme Value Theory of stationary sequences to sampling of large networks. We show that for any general stationary samples (function of node samples) meeting two mixing conditions, the knowledge of bivariate distribution or bivariate copula is sufficient to derive many of its extremal properties. The parameter extremal index (EI) encapsu-Inria