Nicolas Nisse

Ivan Rapaport

Karol Suchan

Distributed computing of efficient routing schemes in generalized chordal graphs *

Keywords: Routing scheme, stretch, chordal graph, distributed algorithm

Efficient algorithms for computing routing tables should take advantage of the particular properties arising in large scale networks. Two of them are of particular interest: low (logarithmic) diameter and high clustering coefficient.

High clustering coefficient implies the existence of few large induced cycles. Considering this fact, we propose here a routing scheme that computes short routes in the class of kchordal graphs, i.e., graphs with no induced cycles of length more than k. In the class of k-chordal graphs, our routing scheme achieves an additive stretch of at most k -1, i.e., for all pairs of nodes, the length of the route never exceeds their distance plus k -1.

In order to compute the routing tables of any n-node graph with diameter D we propose a distributed algorithm which uses messages of size O(log n) and takes O(D) time. The corresponding routing scheme achieves the stretch of k -1 on k-chordal graphs. We then propose a routing scheme that achieves a better additive stretch of 1 in chordal graphs (notice that chordal graphs are 3-chordal graphs). In this case, the distributed computation of the routing tables takes O(min{∆D, n}) time, where ∆ is the maximum degree of the graph. Our routing schemes use addresses of size log n bits and local memory of size 2(d -1) log n bits per node of degree d.

Introduction

In any distributed communication network it is important to deliver messages between pairs of processors. Routing schemes are employed for this purpose. A routing scheme is a distributed algorithm that directs traffic in a network. More precisely, any source node must be able to route messages to any destination node, given the destination's network identifier. When investigating routing schemes, several complexity measures arise. On the one hand, it is desirable to use as short paths as possible for routing messages. Efficiency of a routing scheme is measured in terms of its multiplicative stretch factor (resp., additive stretch factor), i.e., the maximum ratio (resp., difference) between the length of a route computed by the scheme and that of a shortest path connecting the same pair of nodes. On the other hand, as the amount of storage at each processor is limited, the routing information stored in the processors' local memory (the routing tables) must not require too much space with respect to the size of the network. Last but not least, because of the dynamic character of networks, it is important to be able to compute the routing information in an efficient distributed way. While many works propose good tradeoffs between the stretch and the size of routing tables, the algorithms that compute those tables are often impracticable because they are centralized or because of their time-complexity. Indeed, in the context of large scale networks like social networks or Internet, even polynomial time algorithms are inefficient. In this paper, we focus on the tradeoff between the length of the computed routes and the time complexity of computing the routing tables. In the distributed model we consider, the "time" refers to the number of rounds in which O(log n)-bit messages are exchanged between nodes (local computations are "instantaneous").

One way to design efficient algorithms in large scale networks consists in taking advantage of their specific properties. In particular, they are known to have low (logarithmic) diameter and to have high clustering coefficient. Therefore, their chordality (the length of the longest induced cycle) is somehow limited (e.g., see [START_REF] Fraigniaud | Greedy routing in tree-decomposed graphs[END_REF]). That is why, in this paper, we focus on the class of k-chordal graphs. A graph G is called k-chordal if it does not contain induced cycles longer than k. A 3-chordal graph is simply called chordal. This class of graphs received particular interest in the context of compact routing. Dourisboure and Gavoille proposed routing tables of at most log 3 n/ log log n bits per node, computable in time O(m + n log 2 n), that give a routing scheme with additive stretch 2⌊k/2⌋ in the class of k-chordal graphs [START_REF] Dourisboure | Improved compact routing scheme for chordal graphs[END_REF]. Also, Dourisboure proposed routing tables computable in polynomial time, of at most log 2 n bits, but that give an additive stretch k + 1 [START_REF] Dourisboure | Compact routing schemes for generalised chordal graphs[END_REF]. Using a Lexicographic Breadth-First Search (Lex-BFS) ordering (resp., BFS-ordering) of the vertices, Dragan designed a O(n 2 )-time algorithm to approximate the distance between all pairs of nodes up to an additive constant of 1 in n-node chordal graphs, and up to k -1 in, more general, k-chordal graphs [START_REF] Dragan | Estimating all pairs shortest paths in restricted graph families: a unified approach[END_REF]. All these time complexity results consider the centralized model of computation.

In this paper we propose a simpler routing scheme which, in particular, can be quickly computed in a distributed way and achieves good additive stretch for k-chordal graphs. However, the simplicity comes at a price of O(log n) bits per port needed to store the routing tables.

Distributed Model. An interconnection network is modeled by a simple, undirected, connected n-node graph G = (V, E). In the following, D denotes the diameter of G and ∆ denotes its maximum degree. The processors (nodes) are autonomous computing entities with distinct identifiers of size log n bits. We consider an all-port, full-duplex, O(log n)-bit messages, synchronous communication model. That is, links are bidirectional, so, in one communication step, each processor is able to simultaneously send and receive different messages of size O(log n) to and from, respectively, each of its neighbors. By results of Awerbuch and Peleg [START_REF] Awerbuch | Network synchronization with polylogarithmic overhead[END_REF], our algorithms can be adopted to asynchronous networks at a cost of a polylogarithmic overhead. As usual in the context of distributed computation, the time-complexity of algorithms is defined by the required number of communication steps.

Our results. We present a simple routing scheme using a relabeling of the vertices based on particular BFS-trees. Using a Strong-BFS-tree 1 (see Section 2 for the definition), our algorithm achieves an additive stretch k -1 in the class of k-chordal graphs, and using a Maximum Neighborhood BFS-tree (Max-BFS-tree), it achieves an additive stretch 1 in the class of chordal graphs. It uses addresses of size log n bits and local memory of size 2(d -1) log n bits per node of degree d. More precisely, each node must store an interval (two numbers, i.e., 2 log n bits) per port, for all except one port.

The stretches which we achieve are equal to the best ones obtained in previous works. But our algorithm is a (simple) distributed one. It uses O(log n)-bit messages. It computes a relabeling of the vertices and the routing tables in time O(D) when a Strong-BFS-tree is used, and in time O(min{∆D, n}) when a Max-BFS-tree is used. Note that this time-complexity corresponds to the number of communication steps.

In the class of chordal graphs, our results simplify those of Dragan since a Lex-BFS-ordering is more constrained than a Max-BFS-ordering. In particular, the design of a distributed algorithm that computes a Lex-BFS-ordering of the vertices of any n-node graph G in time o(n) is an open problem even if G has small diameter and maximum degree.

Related work. Two kinds of routing schemes have been studied. In the name-independent model, the designer of the routing scheme has no control over the node names (see, e.g., [START_REF] Peleg | A trade-off between space and efficiency for routing tables[END_REF][START_REF] Gavoille | Memory requirements for routing in distributed networks (extended abstract)[END_REF][START_REF] Gavoille | Space-efficiency for routing schemes of stretch factor three[END_REF]). Here we focus on labeled routing, where the designer of the routing scheme is free to name the nodes with labels containing some information about the topology of the network, the location of the nodes in the network, etc. In this context, a routing scheme with multiplicative stretch 4k -5, k ≥ 2, and using Õ(n 1/k ) bits per node 2 in arbitrary graphs is designed in [START_REF] Thorup | Compact routing schemes[END_REF]. In the case of trees, optimal labeled routing schemes using Õ(1) bits per node have been proposed in [START_REF] Fraigniaud | Routing in trees[END_REF][START_REF] Thorup | Compact routing schemes[END_REF]. In [START_REF] Fraigniaud | Routing in trees[END_REF], it is shown that any optimal routing scheme using addresses of log n bits requires Ω( √ n) bits of local memory. Several network classes have been studied, like planar graphs [START_REF] Thorup | Compact oracles for reachability and approximate distances in planar digraphs[END_REF], graphs with bounded doubling dimension [START_REF] Abraham | Routing in networks with low doubling dimension[END_REF], graphs excluding a minor [START_REF] Abraham | Object location using path separators[END_REF], etc. A particular labeled routing scheme is interval routing. Defined in [START_REF] Santoro | Labelling and implicit routing in networks[END_REF], interval routing has received particular interest [START_REF] Gavoille | A survey on interval routing[END_REF]. In such a scheme, the nodes of the network are labeled using integers, and outgoing arcs in a node are labeled with a set of intervals. The set of intervals associated to all the outgoing edges of a node forms a partition of the name range. The routing scheme consists in sending the message through the unique outgoing arc labeled by an interval containing the destination's label. The complexity measure is the maximum number of intervals used in the label of an outgoing arc. An asymptotically tight complexity of n/4 intervals per arc in an n-node network is given in [START_REF] Gavoille | The compactness of interval routing[END_REF]. Moreover, almost all networks support an optimal interval routing scheme using at most 3 intervals per outgoing link [START_REF] Gavoille | The compactness of interval routing for almost all graphs[END_REF]. Specific graph classes have been studied in this context (e.g., bounded treewidth graphs [START_REF] Narayanan | Interval routing on k-trees[END_REF]).

Generalities on BFS-orderings and BFS-trees

In the following, G = (V, E) denotes a connected n-node graph. Let H = (V (H), E(H)) be an induced subgraph of G, i.e., V (H) ⊆ V and

E(H) = {{u, v} ∈ E | u, v ∈ V (H)}. d H (x, y) denotes the distance in H between x, y ∈ V (H). N H (x) denotes the neighborhood of x ∈ V (H) in H. The length |P | of a path P is the number of edges in P . A vertex v ∈ V is simplicial if its neighborhood induces a clique. An ordering {v 1 , • • • , v n } on the vertices of G is called a perfect elimination ordering (PEO) if, for any 1 ≤ i ≤ n, v i is simplicial in G i , where G i is the graph induced by {v i , • • • , v n }.
In the context of a vertex ordering, we denote w < v if w has a smaller index in this ordering. Note that in a PEO, if z < w < v, {z, w} ∈ E and {z, v} ∈ E, then {w, v} ∈ E.

Theorem 1 [FG65] A graph is chordal iff it admits a PEO. 2 The notation Õ() is like O(), up to polylogarithmic factors.

Let r ∈ V . A Breadth-First Search (BFS) ordering of G rooted at r is an ordering of its vertices such that r is the biggest vertex (i.e., the vertex with biggest index in this ordering) and, for any u, v ∈ V (G) \ {r}, v < u (meaning that u has a bigger index than v in the ordering) implies that the biggest neighbor of u is bigger than or equal to the biggest neighbor of v. A Maximum Neighborhood Breadth-First Search (Max-BFS) ordering of G rooted at r is a BFSordering of its vertices with the following additional constraint: for any u, v ∈ V (G) \ {r} with the same biggest neighbor, v < u implies that the number of neighbors of u bigger than u is at least the number of neighbors of v bigger than u. The following theorem will be widely used.

Theorem 2 [BKS05, CK08] A graph G is chordal if and only if any Max-BFS-ordering is a PEO.

Given a BFS-ordering O of the vertices of G, the spanning tree defined by O is the spanning tree obtained by choosing for each vertex, with the exception of the root, its biggest neighbor as the parent. Such a tree defined by a BFS-ordering (resp., by a Max-BFS-ordering) will be called a Strong-BFS-tree (resp., a Max-BFS-tree). Such a tree is rooted at the biggest vertex in the ordering.

Routing Scheme using Strong-BFS and Max-BFS

This section is devoted to presenting a simple routing scheme based on arbitrary spanning trees. We prove that, when Strong-BFS-trees are used, this scheme achieves a good additive stretch in k-chordal graphs, and an improvement of this routing scheme is provided for chordal graphs.

First, let us introduce some notation. Let T be a spanning tree of a graph G. Given x, y ∈ V , T x→y denotes the path in T between x and y. If T is rooted at r ∈ V (T ), its vertices are partitioned into layers: the layer ℓ(v) of a vertex v corresponds to d T (v, r) (which is equal to d G (v, r) if T is a shortest paths tree). If T is rooted, we may say that two vertices are in the same branch if their least common ancestor is equal to one of them. If T is a Strong-BFS-tree and O is an BFS-ordering that defines T , then for given v, w ∈ V , v > w denotes that v has a bigger index than w in the ordering O. Note that, if T is a BFS-tree, then {u, v} ∈ E ⇒ |ℓ(u) -ℓ(v)| ≤ 1. Let us make an observation that will be used in proofs.

Observation 1 Let O be a BFS-ordering that defines a Strong-BFS-tree T rooted at r 0 . Let v ∈ V \ r 0 , l = ℓ(v), and k be any positive integer not greater than l (1 ≤ k ≤ l), and let u be the vertex of maximum index (in O) among the vertices at distance k from v. Then v is in the same branch as u in T .

Proof. Let us proceed by induction on the value of k. If k = 1, then u is the biggest neighbor of v and, by definition of a Strong-BFS-tree, u is the parent of v in T . Now assume that v is in the same branch as u k-1 , the biggest vertex at distance k -1 from v. Since u k-1 has a bigger index than any other vertex b k-1 at distance k -1 from v, its parent u k has an index at least as big as the parent b k of b k-1 .

Finally, given a routing scheme R, Str(R(G), x, y) denotes the difference between the length of the path from x to y computed by R and the distance between x and y in G. The (additive) stretch Str(R(G)) of R in G corresponds to max x,y∈V Str(R(G), x, y).

General Routing Scheme

Let G be a graph and T be any Strong-BFS-tree of G. Roughly speaking, in order to send a message from any source x ∈ V (G) to any destination y ∈ V (G), our routing scheme directs

(R(G 1 , T 1 )) = 3. (b) A chordal graph G 2 and a Max-BFS-tree T 2 such that Str(R(G 2 , T 2 )) = 1. (c) A 2p + 2-chordal graph G 3 and a Max-BFS-tree T 3 such that Str(R(G 3 , T 3 )) = 2p + 1.
the message along the path from x to y in T ; but, if at some step the message could go through an edge e ∈ E(G) \ E(T ) that leads to the branch of T containing y, then it will use such a shortcut. More formally, our routing scheme R(G, T ) is defined as follows.

If x = y, stop. If there is w ∈ N G (x), an ancestor of y in T , choose such a vertex w minimizing d T (w, y); Otherwise, choose the parent of x in T .

For instance, Figure 1 represents three graphs where spanning trees are depicted with bold edges. In Figure 1(a), a message from 1 to 2 will follow the path (1, 4, 6, 7, 5, 2). In Figure 1(b), the same message will follow (1, 4, 5, 2). Let us make some simple remarks.

1. The routing scheme R(G, T ) is well defined. Indeed, the message will eventually reach its destination since its distance to y in T is strictly decreasing at each step. Even if T were not a shortest paths tree, but an arbitrary rooted spanning tree of G, the message would eventually reach its destination. Indeed, in this case, the distance in T between the message and its destination y might not decrease in only one case that may occur at most once: at a step when the message stands at some descendant v of y such that v has a neighbor u (in G) which is an ancestor of y and d T (v, y) ≤ d T (u, y). In that case, the message will be transmitted from v to u, possibly increasing the distance from the message to y in T . It is easy to check that it is the only case when this happens.

Theorem 3 Let k ≥ 3 and let G be a k-chordal graph.

• Let T be any Strong-BFS-tree of G. Then Str(R(G, T )) ≤ k -1.

• Let k = 3 and T be any Max-BFS-tree of G. Then Str(R(G, T )) ≤ 1.

• Both bounds are tight.

Stretch in k-Chordal Graphs

Let k ≥ 3 and let G be a k-chordal graph and T be a (rooted) Strong-BFS-tree of G. Let x, y ∈ V be arbitrary source and destination, respectively. The proof is a case by case analysis to bound Str(R(G, T ), x, y). Let R xy be the route from x to y computed by R(G, T ). In the following, we compare the length of R xy with the length of a shortest path between x and y in G. Several parts of the following discussion are depicted in Figure 2, where bold lines represent edges, thin lines represent paths belonging to T and dotted lines represent paths with edges not necessarily in T .

When the computed route and shortest path are not independent

In this subsection, we prove that it is sufficient to consider x and y with the least common ancestor r 0 such that there is a shortest path P between x and y in G with no internal vertices of P in V (T x→r 0 ) ∪ V (T r 0 →y ), i.e., P and the path with vertex-set V (T x→r 0 ) ∪ V (T r 0 →y ) are independent.

If x is an ancestor or a descendant of y, R xy is the path between x and y in T . Since T is a BFS-tree, this is a shortest path. From now on, we assume that r 0 ∈ V (G), the least common ancestor of x and y, is distinct from x and y. By definition of R(G, T ), R xy either passes through r 0 , or it uses an edge {e, f } ∈ E(G) \ E(T ) with e ∈ V (T x→r 0 ) \ {r 0 } and f ∈ V (T r 0 →y ) \ {r 0 }. I.e., the route R xy from x to y is either T x→e ∪ {e, f } ∪ T f →y , or T x→r 0 ∪ T r 0 →y .

First, we need a technical lemma that shows that the roles of x and y are somehow symmetric.

Lemma 1 Str(R(G, T ), x, y) = Str(R(G, T ), y, x) and the route computed from x to y and the route computed from y to x are almost identical, i.e., they have the same length and differ in at most one vertex.

Proof. Let R xy be the route computed by R(G, T ) from x to y. If R xy passes through r 0 , then there is no edge not in T between a vertex of T x→r 0 and a vertex of T r 0 →y . In this case, R yx also passes through r 0 , and R xy = R yx . Similarly, if R yx passes through r 0 , then R xy = R yx . Now, let us assume that R xy = R yx . Hence, r 0 / ∈ R xy ∪ R yx and both paths contain exactly one edge not in T . This case is illustrated in Figure 2(a). Let {e, f } and {e ′ , f ′ } denote the edges in R x,y \ T and R y,x \ T , respectively. Note that there must be e ′ ∈ V (T e→r 0 ) \ {e} and

f ′ ∈ V (T y→f ) \ {f }. Because T is a BFS-tree, there is d G (r 0 , f ) < d G (r 0 , f ′ ) ≤ d G (r 0 , e ′ ) + 1 ≤ d G (r 0 , e) ≤ d G (r 0 , f ) + 1. Thus, e ′ is the parent of e and f is the parent of f ′ . Therefore, if R xy = R yx , R xy = T x→e ∪ {e, f } ∪ {f, f ′ } ∪ T f ′ →y , and R yx = T y→f ′ ∪ {f ′ , e ′ } ∪ {e ′ , e} ∪ T e→x .
Now, let us have a look at the restrictions we can put on the shortest path that we will compare with the shortest path generated by the routing. Let P 0 be a shortest path in G between x and y. Let y ′ be the first vertex of P 0 in V (T r 0 →y ), and x ′ be the last vertex of P 0 , before y ′ , in V (T x→r 0 ). Let P ′ be the subpath of P 0 between x ′ and y ′ . Notice that P ′ has no internal vertices in V (T x→r 0 ) ∪ V (T r 0 →y ). Because T is a Strong-BFS-tree, P 0 = T x→x ′ ∪ P ′ ∪ T y ′ →y is a shortest path between x and y in G. The following technical lemma restricts our investigation to the case where P 0 = P ′ , i.e., P 0 has no internal vertices in V (T x→r 0 ) ∪ V (T r 0 →y ).

Lemma 2 Either Str(R(G, T ), x, y) = 0, or Str(R(G, T ), x, y) = Str(R(G, T ), x ′ , y ′ ).
Proof. If x ′ = y ′ = r 0 , then it is easy to see that P 0 = R xy . By definition, x ′ and y ′ must be both equal to or different from r 0 . Therefore, let us assume both are different from r 0 . Recall that R xy = T x→e ∪ {e, f }∪ T f →y , or R xy = T x→r 0 ∪ T r 0 →y . In the second case, we set e = f = r 0 . The proof is a case by case analysis according to the relative positions of x ′ and e in T x→r 0 , and of y ′ and f in T r 0 →y .

If

T x→e ⊆ T x→x ′ and T f →y ⊆ T y ′ →y , then Str(R(G, T ), x, y) = 0 because |P ′ | ≥ 1. Let us assume that T x→x ′ ⊂ T x→e and T f →y ⊂ T y ′ →y . This case is illustrated in Figure 2(b). Note that in this case e = f , therefore {e, f } ∈ E(G). Let a = |T x→e | -|T x→x ′ | > 0, and let b = |T y ′ →y | -|T f →y | > 0.
We study the layers of x, x ′ , y and y ′ to prove that Str(R(G, T ), x, y) = 0. Let L = ℓ(e) be the layer of e. Then, ℓ(x

′ ) = L + a. Because {e, f } ∈ E(G) and T is a BFS-tree, L -1 ≤ ℓ(f ) ≤ L + 1. Therefore, L -1 -b ≤ ℓ(y ′ ) ≤ L + 1 -b. However, because T is a BFS-tree, we must have ℓ(x ′ ) ≤ ℓ(y ′ ) + |P ′ |. Thus, ℓ(x ′ ) ≤ L + 1 -b + |P ′ |. Finally, we get that a + b -1 ≤ |P ′ |. Since b > 0, then a ≤ |P ′ |. To conclude, let us observe that |R xy | = |T x→x ′ | + a + 1 + |T y ′ →y | -b = |P | -|P ′ | + a + 1 -b ≤ |P |. Hence, Str(R(G, T ), x, y) = 0.
If T y→y ′ ⊂ T f →y and T e→x ⊂ T x ′ →x , the reasoning given above shows that Str(R(G, T ), y, x) = 0. By Lemma 1, we get that Str(R(G, T ), x, y) = 0.

Finally, if

T x→x ′ ⊆ T x→e and T y ′ →y ⊆ T f →y , the route computed by R(G, T ) from x ′ to y ′ is clearly T x ′ →e ∪ {e, f } ∪ T f →y ′ . Moreover, Str(R(G, T ), x, y) = |T x ′ →e ∪ {e, f } ∪ T f →y ′ | -|P ′ | = Str(R(G, T ), x ′ , y ′ ).

When the computed route and shortest path are independent

From Lemma 2, it remains to consider the case where there exists an xy-shortest path P in G with no internal vertices of P in V (T r 0 →y ) ∪ V (T x→r 0 ), where r 0 is the least common ancestor of x and y in T (cf. Figures 2(c)-2(f)). Basically, the proof proceeds by analyzing the distances in the cycle T x→r 0 ∪ T r 0 →y ∪ P and finding convenient chords in it.

Let S ⊆ V , N G (S) denotes the set of all vertices in V \ S with at least one neighbor in S.

Claim 1 If r 0 / ∈ N G (P ), there exist u in T x→r 0 , u = x, and v in T y→r 0 , v = y, such that u, v ∈ N G (P ). Moreover, if G is chordal, u and v may be chosen adjacent: {u, v} ∈ E(G)\E(T ).

Proof. Since r 0 / ∈ N G (P ), the parents in T of x and y satisfy the first part of the claim. For the second part, consider that G is chordal. Let C r 0 be the connected component of G \ N G (P ) that contains r 0 . Let S = N G (C r 0 ) and let C P be the connected component of G \ S that contains P . Notice that S ⊆ N G (P ) and it is an inclusion-minimal xr 0 -separator.

We can pick vertices u and v from the intersections of S with T x→r 0 and T r 0 →y , respectively (see Fig. 2(c)). Since S is a minimal separator in a chordal graph, S induces a clique [START_REF] Golumbic | Algorithmic Graph Theory and Perfect Graphs[END_REF], Thus {u, v} ∈ E(G). Finally, {u, v} / ∈ E(T ) because the opposite would imply that u or v is the a common ancestor of x and y, i.e., r 0 ∈ {u, v}, a contradiction since u, v ∈ S. ⋄ Lemma 3 Let k ≥ 3. Let G be a k-chordal graph and let T be a Strong-BFS-tree of G. Then, for any x, y ∈ V (G), Str(R(G, T ), x, y) ≤ k -1.

Proof. Recall that we only need to consider the following case: r 0 is the least common ancestor of x and y, and there exists a shortest xy-path P with no internal vertex in T x→r 0 ∪ T r 0 →y . If the route R xy computed by R(G, T ) takes a shortcut, this edge is denoted {e, f }.

There are two cases to be considered. We first assume that r 0 is not in the neighborhood of P , N G (P ) (cf., Figures 2(d Let us choose u and v as defined in Claim 1 and such that d G (u, r 0 ) + d G (r 0 , v) is minimum. We can assume that e ∈ T r 0 →u and f ∈ T v→r 0 (possibly e = f = r 0 ). Indeed, we can show the following.

Claim 2 Let {e, f } be the "shortcut" used in R xy . If T r 0 →e or T r 0 →f has an internal vertex a in N G (P ), then Str(R(G, T ), x, y) ≤ 2.

Proof. Indeed, suppose that the path T r 0 →e has an internal vertex a in N G (P ). In other words, there is a vertex a in the neighborhood of P , in the same branch as e and d G (x, a) > d G (x, e). Let z denote a neighbor of a in P . Since T is a BFS-tree,

d G (x, z)+1 ≥ d G (x, a), d G (r 0 , f )+1 ≥ d G (r 0 , e), and d G (r 0 , a) + 1 + d G (z, y) ≥ d G (r 0 , f ) + d G (f, y). So d G (x, z) ≥ d G (x, e), d G (r 0 , f ) ≥ d G (r 0 , a), d G (r 0 , f ) + 1 + d G (z, y) ≥ d G (r 0 , f ) + d G (f, y), and 1 + d G (z, y) ≥ d G (f, y). Finally, d G (x, z) + 1 + d G (z, y) + 1 ≥ d G (x, e) + d G (f, y) + 1. Thus Str(R(G, T ), x, y) ≤ 2. By a similar argument, b, an internal vertex of T r 0 →f in N G (P ) results in Str(R(G, T ), x, y) ≤ 2.
⋄ Thus, it is left to consider the case where P has no such neighbors, and e ∈ T r 0 →u and f ∈ T v→r 0 . In particular, this means that there are no edges between a vertex in T e→u and

T v→f but {e, f }. Let u ′ ∈ N G (u) ∩ P and v ′ ∈ N G (v) ∩ P such that d G (u ′ , v ′ ) is minimum.
In the following we assume that u ′ is between x and v ′ in P and let P = P 1 ∪ P 2 ∪ P 3 , where P 1 is the subpath of P between x and u ′ , P 2 is the subpath of P between u ′ and v ′ , and P 3 is the subpath of P between v ′ and y. Otherwise, the proof is similar by setting P 1 to be the subpath of P between x and v ′ , P 2 is the subpath of P between v ′ and u ′ , and P 3 is the subpath of P between u ′ and y.

By the choice of u, u ′ , v, v ′ , the cycle 2(d)). Therefore, let us consider the case when |P 2 | = 0, i.e., u ′ = v ′ (cf., Figure 2(e)). We first consider the case when u > v (in the BFS-ordering defining T ). There must be |T v→y | ≤ |P 3 |. Indeed, suppose |T v→y | = |P 3 | + 1. This means that v is not the biggest vertex at distance d G (v, y) from y (u is bigger) and, by Observation 1, y cannot be in the same branch as v in T -a contradiction. Hence,

C = {u, u ′ } ∪ P 2 ∪ {v ′ , v} ∪ T v→f ∪ {f, e} ∪ T e→u has no chord. Thus, |C| = 3 + |P 2 | + |T v→f | + |T e→u | ≤ k. Because T is a BFS-tree, |T x→u | ≤ 1 + |P 1 | and |T v→y | ≤ 1 + |P 3 |. It follows that |R xy | = |T x→u | + |T v→y | + 1 + |T v→f | + |T e→u | ≤ |P 1 | + |P 3 | + k -|P 2 |. If |P 2 | > 0, |R xy |-|P | = Str(R(G, T ), x, y) ≤ k -2|P 2 | < k -1 (cf., Figure
|R xy | = |T x→u |+|T v→y |+1+|T v→f |+|T e→u | ≤ |P 1 |+|P 3 |-1+k-|P 2 | = k+|P 1 |+|P 3 |-1 = |P |+k-1.
Now, let us consider the case when u < v. Exchanging the roles of x and y, by symmetry, the above reasoning shows that Str(R(G, T ), y, x) ≤ k -1. Now, applying Lemma 1, we again get that Str(R(G, T ), x, y) ≤ k -1.

To conclude, let us assume that r 0 ∈ N G (P ) (cf., Figure 2(f)). Let z ∈ N G (r 0 ) ∩ P . Let P = P 1 ∪ P 2 where P 1 is the subpath of P between x and z, and P 2 is the subpath of P between z and y. Because T is a Strong-BFS-tree,

|T x→r 0 | ≤ 1 + |P 1 | and |T r 0 →y | ≤ 1 + |P 2 |. Therefore, |R xy | ≤ |T x→r 0 | + |T r 0 →y | ≤ |P 1 | + |P 2 | + 2 ≤ |P | + k -1 (because k ≥ 3).
It is important to note that the previous result is valid for any Strong-BFS-tree. However, it is easy to observe that the inequality given by Lemma 3 is optimal. Indeed, Figure 1(c) represents a k-chordal graph with k = 2p + 2 (p ≥ 1) and a Strong-BFS-tree T (that actually is a Max-BFS-tree) such that Str(R(G, T ))) = 2p + 1 = k -1: a message from 1 to 2 will pass through the edge {3p + 3, 3p + 4}.

Lemma 3 gives that, for any chordal graph G and a Strong-BFS-tree T , Str(R(G, T )) ≤ 2. The following lemma states that in the class of chordal graphs we can improve the stretch down to 1 by using a "better" BFS-tree, i.e., a Max-BFS-tree.

Lemma 4 Let G be a chordal graph and let T be a spanning tree defined by any Max-BFSordering. Then, for any x, y ∈ V (G), Str(R(G, T ), x, y) ≤ 1.

Proof. By Lemma 2, it is enough to consider the case when r 0 is the least common ancestor of x and y, P is a shortest path between x and y that is independent from T x→r 0 ∪ T r 0 →y , i.e., the only common vertices of P and T x→r 0 ∪T r 0 →y are x and y. There are two cases to be considered.

1. We first assume that r 0 / ∈ N G (P ).

The proof consists of several simple claims. This case is illustrated in Figures 2(c Let us choose u and v as defined in Claim 1. Note that {u, v} ∈ E because G is chordal. Therefore, by definition of the routing scheme, and because T is a shortest paths tree:

Claim 3 |R xy | ≤ |T x→u | + 1 + |T v→y |.
Claim 4 u and v have a common neighbor z in P .

Proof. Let u ′ and v ′ be neighbors of u and v respectively in P . u ′ and v ′ are chosen in order to minimize the length of the subpath P ′ of P between u ′ and v ′ . Consider the cycle {u, v} ∪ {v, v ′ } ∪ P ′ ∪ {u ′ , u}. Note that it is an induced cycle and, because G is chordal, has to be a triangle. Thus u ′ = v ′ = z. ⋄ Let z be a common neighbour of u and v in P .

Claim 5 If z ∈ {x, y}, Stretch(R(G, T ), xy) ≤ 1.
Proof. We assume z = x (the case z = y is symmetric). This case is illustrated in Figure 2

Because {u, v} ∈ E, ℓ(u) -1 ≤ ℓ(v) ≤ ℓ(u) + 1. If ℓ(v) = ℓ(u) + 1, then ℓ(y) = ℓ(v) + |T v→y | = ℓ(u) + 2 + |P |. Besides, ℓ(y) ≤ ℓ(x) + |P |. Since ℓ(x) = ℓ(u)
+ 1, we get a contradiction. Now, let us assume ℓ(v) = ℓ(u) -1. In this case, v > u and v should be the parent of x in a Strong-BFS-tree T , but we have {u, x} ∈ E(T ) and {v, x}

∈ E \ E(T ) -a contradiction. Finally, if ℓ(v) = ℓ(u), let w ′ be a vertex in P \ {x, y}. If ℓ(w ′ ) ≤ ℓ(x), then ℓ(v)+|T v→y | = ℓ(v)+1+|P | = ℓ(y) ≤ ℓ(w ′ )+|P |-1 ≤ ℓ(x)+|P |-1 = ℓ(u)+|P | = ℓ(v)+|P |
-again a contradiction. Therefore, ℓ(w ′ ) > ℓ(x) for any w ′ ∈ P \ {x, y}. This implies that, for any w ′ ∈ P \ {x, y}, ℓ(w ′ ) -ℓ(v) > 1 and thus, {v, w ′ } / ∈ E (property of levels in a BF S tree). Also, for any vertex w ′′ ∈ T v→y \ {w, v}, {v, w ′′ } / ∈ E because T v→y is a shortest path. To conclude, let us consider the cycle P ∪ T y→v ∪ {v, x}. Since G is chordal, we must have {x, w} ∈ E. ⋄ Therefore, let us assume that z / ∈ {x, y}.

Moreover, let us assume that u > v (otherwise, we prove symmetrically that Str(R(G, T ), y, x) ≤ 1 and Lemma 1 allows to conclude).

Let P 1 be the subpath of P between x and z, and let P 2 be the subpath of P between z and y. Because T x→u and T v→y are shortest paths, we get: Proof. As above, z > w * would contradict x being in the same branch as w * . To conclude, if {u, z} / ∈ E(T ), then the parent of z in T must be strictly bigger than u, which implies that z > w * , a contradiction.

Claim 6 |T x→u | ≤ |P 1 | + 1 and |T v→y | ≤ |P 2 | + 1. Claim 7 |T v→y | ≤ |P 2 |.
⋄ Claim 9 If |T x→u | = |P 1 | + 1,
| = |P 1 | + 1, then v > w * . Proof. Let v ′ be the parent of v in T . We first prove that v ′ > u. Since u > v and {u, v} ∈ E \ E(T ), ℓ(u) ≤ ℓ(v) ≤ ℓ(u) + 1. If ℓ(u) = ℓ(v), then ℓ(u) = ℓ(v ′ ) + 1 and v ′ > u.
Otherwise, ℓ(u) = ℓ(v ′ ), u and v ′ are both adjacent to v and {v ′ , v} ∈ E(T ). Hence, v ′ > u. Since the parent of v is bigger than the parent of w * , we obtain that v > w * . ⋄

Claim 11 If |T x→u | = |P 1 | + 1, then {w * , v} ∈ E.

2.

Let us now assume that r 0 ∈ N G (P ).

Let z ∈ N G (r 0 ) ∩ P , let u be the child of r 0 in T r 0 →x and let v be the child of r 0 in T r 0 →y . This case is illustrated in Figure 2(h). Therefore, let us assume z / ∈ {x, y}. Let P 1 be the subpath of P between x and z, and let P 2 be the subpath of P between z and y. Because T x→r 0 and T r 0 →y are shortest paths: It is easy to observe that the inequality given by Lemma 4 is optimal. Indeed, consider Figure 1(b) and the route between 1 and 2. The above discussion and lemmata prove Theorem 3.

Claim 12 |T x→r 0 | ≤ |P 1 | + 1 and |T r 0 →y | ≤ |P 2 | + 1. Because |R xy | ≤ |T x→r 0 | + |T

Distributed algorithm

In this section we present a simple distributed algorithm that computes the routing tables sufficient for the execution of the routing scheme described in the previous section. Recall that the time-complexity of the algorithm is defined by the number of communication steps. We neglect the time needed for computations performed locally, since passing messages between the nodes is considered several orders of magnitude slower. To reinforce this point of view let us point out that the local computations used by the algorithm are not complex. Actually, the most costly computation performed by the nodes during each step of the execution of our distributed algorithm consists in lexicographically ordering ∆ labels of at most D integers of log D bits which can be done in time ∆D.

Before describing the algorithms for computing routing schemes, let us give an informal description of the procedure. The computation of routing schemes consists of three phases. The first two of them aim at building a Strong/Max-BFS-tree T . For this purpose, the first phase chooses an arbitrary vertex as the root and gives to any vertex its layer, i.e., its distance to the root. The second phase is divided into rounds. During the i th round, each vertex in layer i decides its parent in layer i -1 and the vertices in layer i -1 order their children. The ordering process is done simultaneaously for all vertices in layer i -1, in time O(1) in the case of a Strong-BFS-tree, resp., in time O(∆) in the case of a Max-BFS-tree. Then, during the third phase, each vertex x is assigned an integral label P (x) that corresponds to its position in a DFS postorder traversal of T . It is easy to check that it gives a Strong/Max-BFS-ordering of G. Moreover, x learns I(x), the interval corresponding to the labels of its descendants in T (including x). P (x) is used as the identifier of x in the routing scheme. At each vertex y, for every neighbor x of y (except for the parent of y), the edge yx is labeled with I(x).

Let us describe the algorithm in detail.

1 st Phase. The first phase chooses an arbitrary vertex r ∈ V (G) as the root and gives to each vertex its layer, i.e., its distance from r. Moreover, each vertex informs its neighbors of its own layer. This trivially takes at most D + 1 steps by broadcasting a counter initially set to 1 by the root. Now, if each vertex chooses an arbitrary neighbor in the lower layer as the parent, the obtained graph is a BFS-tree. However, when Strong-BFS-trees or Max-BFS-tree are considered, a particular neighbor in the lower layer must be chosen as parent.

2 nd Phase. The second phase aims at determining an appropriate parent for each vertex. For this purpose, we assign an ordering on the vertices based on the following labeling: the root receives an empty label and any vertex v ∈ V (G) in the layer i ≥ 1 will eventually have a full label label(v), where label(v) is a sequence of i integers that consists of the full label of its parent u concatenated with the integer p that indicates that v is the p th child of u. In the following, a partial label is a prefix of a full label that will eventually be computed. Note that, if two nodes u and v have same partial label with size i ≥ 1, i.e, the first i integers of their full labels are equal, then u and v have the same ancestors in layers 0 to i. Roughly, the label of a vertex v defines a path from the root to v. The labels will be constructed gradually, in a way that each vertex will be aware of the current (partial) labels of its neighbors. Notice that the lexicographic ordering of full labels gives the inverse of the Strong-BFS-ordering (or Max-BFS-ordering) under construction. Transforming it into integer numbers ranging from n down to 1 can be easily computed once we have fixed T and ordered the children of each node (see the third phase).

The second phase is devided into D rounds. At the end of the i th round, vertices of layer i ≥ 1 have learnt their full label, while the vertices of layer i + 1 know their partial label of length i. The vertices in layer i + 1 also know the full label of their neighbors in layer i and they can choose as parents the one among these neighbors with smallest full label. Then, during the (i + 1) th round, the vertices of layer i order their children so that any vertex v in layer (i + 1) learns the last integer of its (full) label. This takes time O(1) when Strong-BFS-trees are considered, and time O(∆) when Max-BFS-trees are considered. Once, this integer has been computed, the vertex v propagates it toward layers > i + 1. This spreading of messages of size O(log n) bits is done simultaneously to next rounds, in such a way that when round j occurs, the vertices in layer j actually know their partial label of length j -1.

Spreading of labels. Let us describe how to spread the labels of the vertices efficiently, that is, using messages of size O(log n) bits and in time O(D).

It is easy to check that the following process ensures that k ≥ 0 steps after the end of the i th round (i ≥ 1), the vertices in layer i + k + 1 have learnt their partial label of length i and the partial label of length i of their neighbors in layers i + k + 1 and i + k. In particular, this ensures that before the beginning of roud i + 1, all vertices of layer i + 1 know the full label of their neighbors in layer i and their own partial label of length i.

For any i ≥ 1, each vertex v in the layer i maintains a subset P P (v) (for potential parent) of its neighbors in layer i -1 that initially contains all these neighbors. Roughly, at any step of the algorithm, all vertices in P P (v) will have the same partial label which correponds to the partial label of v at this step. When all vertices in P P (v) receive a new integer to be added in their partial label, they forward it to v that only keeps as potential parents the vertices in P P (v) with smallest partial label.

More formally, once a vertex v in layer i has received its full label with last integer p v (that corresponds to its position among its siblings), it transmits the pair (i, p v ) to all its neighbors. Then, its neighbors in layer i + 1 must propagate this information toward the layers j > i + 1. To avoid flooding, we proceed as follows, any vertex considering only the messages coming from its potential parents. When a vertex v in the layer j > i has received such a message (i, p u ) from all of its neighbors u in P P (v), this means that, at this step, the partial labels of the vertices in P P (v) have length i and differ only in their i th element. That is to say, there is a sequence of i -1 integers L such that for any u ∈ P P (v), the label of u is L concatenated with p u . Moreover, L is the partial label of v at this step. Then, v keeps in P P (v) the vertices that have the smallest partial label (in the lexicographic ordering): let p be smallest integer such that v has received (i, p) from its neighbors in P P (v), then v keeps in P P (v) all the vertices with new partial label L concatenated with p. Then, v adds p in its own partial label and transmits the corresponding pair (i, p) to all its neighbors. Moreover, receiving such a message from any neighbor u ′ , v adds p to the locally stored (partial) label of u ′ . Proceeding in this way, once all vertices v in layer i have received a label, every vertex in layer i + 1 knows its potential parents, i.e., its neighbors in P P (v), and the corresponding label.

Ordering of children in Strong-BFS-tree. Once each vertex in layer i has chosen a parent, the vertices in layer i -1 establish an ordering on their children. If we want to obtain a Strong-BFS-tree without additional properties, any ordering is valid. Therefore, each vertex in the layer i -1 arbitrarily orders its children and sends them their position in this ordering. Then, each vertex in layer i has a full label. This takes one step per layer, i.e., time O(D) in total.

Ordering of children in Max-BFS-tree. In this case, each vertex in layer i -1 will order its children according to the number of neighbors with smaller labels they have. In other words, each vertex in layer i -1 orders its children according to the number of their neighbors that will have larger numbers in the final ordering. Notice that as soon as the vertices of layer i have chosen their parent and broadcasted them to their neighbors, a vertex v in layer i only needs to learn its position in the ordering relatively to its siblings in T .

Note that, by definition of Max-BFS, the order of the children of a vertex in layer i -1 does not depend on the order of the children of another vertex in layer i -1. The children of vertices in layer i -1 appear as consecutive blocks in layer i and it is only important to know which blocks appear earlier and which ones appear later in the ordering. Therefore, the children of each vertex u from layer i -1 can be ordered in parallel, based on the number of their neighbors in layer i -1, or in layer i and with a father placed before u in layer i -1.

A vertex u in layer i -1 orders its children as follows. Let us assume u has already ordered its first p children (p ≥ 0). These neighbors of u have full labels while remaining neighbors of u only have partial label. By definition of a Max-BFS-ordering, u chooses its p + 1 st child v as the one with the biggest number of neighbors with bigger labels (than v), i.e., v is the child with the biggest number of neighbors that either have a parent bigger than u or that are siblings of v with a full label. v receives p + 1 from u and completes its label (that becomes full). Then, v informs its siblings that it has received a full label, and each child of u updates the number of its neighbors that already have full labels. In this way, u orders its d(u) children in O(d(u)) time. So, in total, this step is executed in O(∆) time per layer. Therefore, in at most O(min{∆D, n}) steps, every vertex has chosen a unique vertex as its parent and the tree T rooted in r is well defined.

3 rd Phase. The third phase consists in assigning to each vertex v his position in the ordering and the interval of positions of vertices that belong to T v , the subtree of T rooted in v. It is easy to do so by two stages, the first one consisting in propagation of messages from the leaves toward the root and the second one from the root toward the leaves. During the first stage, the leaves of T send 1 to their parents, and every vertex u with children v 1 , • • • , v r receives from v i (1 ≤ i ≤ r) the number ℓ i of vertices belonging to the subtree of T rooted in v i and sends to its parent 1 + i≤r ℓ i . During the second stage, the root is assigned the position n and the interval [1, . . . , |V (G)|]; each vertex v takes the last position in the interval it has received and partitions the rest into subintervals corresponding to each of its children. It is easy to check that the resulting ordering corresponds to a DFS postorder traversal of T . This phase takes at most 2D steps. The discussion of this section can be summarized with the following theorem.

Theorem 4 The distributed protocol described above computes routing tables of O(∆ log n) bits per node for any n-node network G with diameter D and maximum degree ∆. Moreover, our distributed protocol is executed in time O(D) if the tree T (implicit in the table) is an arbitrary Strong-BFS-tree, and in time O(min{∆D, n}) if T is a Max-BFS-tree. The routing tables produced by our prococol are the one used by the routing scheme described in Section 3.

The distributed protocol described in this section has been implemented on a simulator and has been tested for various topologies [START_REF] Hogie | Alusim: simulating routing schemes on large-scale topologies[END_REF]. For each topology, and for each networks' size in {500, 1000, • • • , 3500}, 100 random instances of networks with the considered topology have been tested. The results obtained via the simulations correspond to theoretical results. In particular, an empiric argument for the correctness of the implementation is that, in chordal graphs, the additive stretch never exceeds 1. Non surprisingly, in grid topologies which include many big holes, our protocol achieves poor performances (roughly, an average additive stretch of n is obtained for the n × n grids). Further, the protocol behaves very well on power law networks (generated using the preferential attachment method) which are considered to be good representations of real world complex networks [START_REF] Barabasi | Emergence of scaling in random networks[END_REF]. Indeed, in this kind of networks, our protocol achieves an average (resp., maximum) additive stretch less than 1 (resp., 5).

Figure 1 :

 1 Figure 1: Lemmata 3 and 4 give optimal bounds. (a) A 4-chordal graph G 1 and a Strong-BFS-tree T 1 such that Str(R(G 1 , T 1 )) = 3. (b) A chordal graph G 2 and a Max-BFS-tree T 2 such that Str(R(G 2 , T 2 )) = 1. (c) A 2p + 2-chordal graph G 3 and a Max-BFS-tree T 3 such that Str(R(G 3 , T 3 )) = 2p + 1.

  ), 2(e)).

Figure 2 :

 2 Figure 2: Illustrations of different cases when bounding Str(R(G, T ), x, y). 2(a) Lemma 1, R xy = R yx . 2(b) Lemma 2, T x→x ′ ⊂ T x→e and T f →y ⊂ T y ′ →y . 2(c) Claims 1, 8, 9 and 10. 2(f) Lemma 3, r 0 ∈ N G (P ). 2(d) Lemma 3, r 0 / ∈ N G (P ), first item, case |P 2 | > 0. 2(e) Lemma 3, r 0 / ∈ N G (P ), first item, case |P 2 | = 0. 2(g) Claim 5. 2(h) Lemma 4, r 0 ∈ N G (P ).

Proof.

  Indeed, suppose |T v→y | = |P 2 | + 1. This means that v is not the biggest vertex at distance d G (v, y) from y (u is bigger) and, by Observation 1, y cannot be in the same branch as v in T -a contradiction. ⋄ Let w * be the child of u in T u→x . Next three claims are illustrated in Figure 2(c). Claim 8 If |T x→u | = |P 1 | + 1, then w * > z and {u, z} ∈ E(T ).

Proof.

  By Claims 8 and 10, we get v > w * > z. Moreover, by Claim 9, {w * , z} ∈ E. By definition of z, {z, v} ∈ E. By Theorem 2, the ordering defined by a Max-BFS is a PEO, therefore, {w * , v} ∈ E. ⋄ We are now ready to prove first case of Lemma 4. If the common neighbor z of u and v in P (z exists by Claim 4) is either x or y, then the result follows from Claim 5. Otherwise, by Claim 6, |T x→u | ≤ |P 1 | + 1. By Claims 7 and 11, it follows that |T v→y | ≤ |P 2 | and, either |T x→u | ≤ |P 1 |, or |T x→u | = |P 1 | + 1 and {w * , v} ∈ E. If |T x→u | ≤ |P 1 |, by Claim 3, Str(R(G, T ), x, y) ≤ 1. Otherwise, because {w * , v} ∈ E and by definition of the routing scheme, |R xy | ≤ |T x→w * | + 1 + |T v→y | and then |R xy | ≤ |P 1 | + 1 + |P 2 |, i.e., Str(R(G, T ), x, y) ≤ 1. This proves the first case of Lemma 4.

First, let us

  consider the case z = x (or symmetricaly z = y). Because T r 0 →y is a shortest path, |T r 0 →y | ≤ 1 + |P |. Moreover, |R xy | ≤ 1 + |T r 0 →y |. Hence, if |T r 0 →y | ≤ |P |, Str(R(G, T ), x, y) ≤ 1. Therefore, let us consider the case when |T r 0 →y | = 1 + |P |. In particular, y = v because |P | > 0. If {z, v} ∈ G, |R xy | ≤ 1 + |T v→y | ≤ 1 + |P |, and Str(R(G, T ), x, y) ≤ 1. It remains to prove that {z, v} ∈ G. Let us consider the cycle C = {r 0 , v}∪T v→y ∪P ∪{z, r 0 }. For any vertex t ∈ C \{r 0 , z, v}, because |T r 0 →x | = |P |+1, ℓ(t) ≥ ℓ(r 0 )+2. Therefore {v, t} / ∈ E. By chordality of G and because |C| ≥ 4, {z, v} ∈ E.

  r 0 →y | and by Claim 12: Claim 13 If |T x→r 0 | ≤ |P 1 | or |T r 0 →y | ≤ |P 2 |, Str(R(G, T ), x, y) ≤ 1. From now on, let us assume that |T x→r 0 | = |P 1 | + 1 and |T r 0 →y | = |P 2 | + 1. Similarly to the above case when z = x, we prove that {z, v} ∈ E and {z, u} ∈ E. By mimicking the proof of Claim 8, we get that z < min{u, v}. By Theorem 2, the ordering defined by a Max-BFS is a PEO, therefore, {u, v} ∈ E. Then, |R xy | ≤ |T x→u | + 1 + |T v→y |, and Str(R(G, T ), x, y) ≤ 1. This concludes the proof of the second case of Lemma 4.

  (g). By construction, there is |R xy | ≤ 1 + |T v→y | and |T v→y | ≤ 1 + |P |. If |T

v→y | ≤ |P |, Str(R(G, T ), x, y) ≤ 1. So we only have to consider the case |T v→y | = |P |+1. Let w be the child of v in T v→y . If {x, w} ∈ E, then |R xy | ≤ 1 + |T w→y | ≤ |P | + 1 and Str(R(G, T ), x, y) ≤ 1. We prove that {x, w} ∈ E.

  then {w * , z} ∈ E. Proof. Let us consider the cycle C = {u, w * } ∪ T w * →x ∪ P 1 ∪ {z, u}. For any vertex t ∈ C \ {u, w * , z}, because |T x→u | = |P 1 | + 1, ℓ(t) ≥ ℓ(u) + 2. Therefore {u, t} / ∈ E. By chordality of G and because |C| ≥ 4, {w * , z} ∈ E. ⋄ Claim 10 If |T x→u

The name BFS-tree is often used in distributed computing literature to denote any shortest paths tree. To emphasize the particular properties of BFS-trees that are used in this work, the authors decided to add the prefix "Strong" in the name, even though the BFS-trees found in textbooks on algorithms usually are Strong-BFS-trees in this sense.

Once a spanning tree T rooted at an r ∈ V (G) has been defined, this scheme can be efficiently implemented. It is sufficient to label the vertices in a way where any rooted subtree of T corresponds to a single interval. Each vertex u stores the interval corresponding to the subtree of T rooted at v, for any neighbor v of u but its parent. Then, the routing function chooses the port corresponding to the inclusion-minimal interval containing the destination's address, and it chooses the parent of the current location if no such interval exists. Note that this is not an interval routing scheme because some intervals may be contained in others.

3. Since we assume that T is a BFS-tree, it is easy to see that the route computed by the routing scheme R(G, T ) between two arbitrary nodes contains at most one edge that is not an edge of T . Indeed, after having taken such a shortcut, the message reaches y by following the path in T , which is a shortest path in G since T is a BFS-tree. Subsection 3.2 below is devoted to proving the following theorem.

* Partially supported by programs of Conicyt Chile: Anillo ACT88 (K.S.), Basal-CMM (I.R. and K.S.), Fondecyt 1090156 (I.R.), Fondecyt 11090390 (K.S.), Instituto Milenio ICDB (I.R.) and the European Commission through the EULER project (Grant No.258307) part of the Future Internet Research and Experimentation (FIRE) objective of the Seventh Framework Programme (FP7). (N.N.).

Open Problems

Many questions remain open in this study. In particular: Is it possible to design a routing scheme achieving the same stretch and time-complexity but using smaller routing tables? What can we say about the stretch if few large cycles are allowed? Routing schemes in dynamic networks, i.e., when nodes are free to leave or to arrive in the network at any time, are needed. Fault-tolerant and self stabilizing algorithms to compute routing tables would be of much interest.

The authors would like to thank Feodor F. Dragan for valuable discussions during his stay at CMM, Universidad de Chile.


Lélia Blin

Janna Burman

Nicolas Nisse

Brief Announcement: Distributed Exclusive and Perpetual Tree Searching ⋆

We tackle a version of the well known graph searching problem where a team of robots aims at capturing an intruder in a graph. The robots and the intruder move between the graph nodes. The intruder is invisible, arbitrary fast, and omniscient. It is caught whenever it stands on a node occupied by a robot, and cannot escape to a neighboring node. We study graph searching in the CORDA model of mobile computing: robots are asynchronous and perform cycles of Look-Compute-Move actions. Moreover, motivated by physical constraints and similarly to some previous works, we assume the exclusivity property, stating that no two or more robots can occupy the same node at the same time. In addition, we assume that the network and the robots are anonymous. Finally, robots are oblivious, i.e., each robot performs its move actions based only on its current "vision" of the positions of the other robots. Our objective is to characterize, for a graph G, a set of integers such that for every integer k in the set, perpetual graph searching can be achieved by a team of k robots starting from any k distinct nodes in G. One of our main results is a full characterization of this set, for any asymmetric tree. Towards providing a characterization for all trees, including trees with non-trivial automorphisms, we have also provided a set of positive and negative results, including a full characterization for any line. All our positive results are based on the design of graph searching algorithms.

Introduction

This BA announces a work that aims at understanding the algorithmic power and limitation of computing with autonomous mobile robots. The literature dealing with this objective has considered different kinds of coordination tasks, including pattern formation, robot gathering, and graph exploration. Each of these tasks involves complex coordination protocols for the robots, whose complexity depends on the capabilities of the robots in terms of perception of their environment, individual computational power, and communication. The CORDA model has been introduced for capturing the essence of mobile computing. It focusses on the asynchronous nature of the actions performed by the robots, and on the limitations caused by the absence of communication between them.

In this model, robots are endowed with visibility sensors allowing each robot to perceive the positions of all the other robots. They operate in asynchronous Look-Compute-Move action cycles. During its look action, a robot perceives the relative positions of the other robots; during the compute action, it executes some individual deterministic computation whose input is the set of the latest perceived positions of the other robots (which may have changed since they were measured); finally, during the move action, the robot changes its position according to its computation.

The coordination of autonomous deterministic mobile robots has been first studied in continuous environments (e.g., the 2-dimensional Euclidean space).

In the discrete CORDA model, it is not clear whether even simple coordination tasks can be achieved. As a consequence, most of the literature on robot computing assume additional hypotheses, including presence of identities enabling to distinguish robots, capacity to store the sequence of all previous perceived positions and moves, sense of direction, and the ability to construct towers of robots (i.e., to bring several robots at a node). We have proved that a very complex task, like graph searching, can be solved without any of these assumptions. We assume that the robots have no identities, they are memoryless (i.e., stateless), have no sense of direction, and two or more robots cannot occupy the same position at the same time (exclusivity property). We call min-CORDA this essential CORDA model.

Contributions

For any graph G, let fs(G) denote the set of robot team sizes such that k ∈ fs(G) if and only if distributed graph searching in G can be achieved by a team of k robots starting from any k (distinct) nodes. Our main result consists in a characterization of fs(T ), for any asymmetric tree T . Let xs(T ) denote the exclusive search number, i.e., the minimum number of robots to be used for capturing the intruder in a centralized setting satisfying the exclusivity property. By definition, k ∈ fs(T ) implies that k ≥ xs(T ). We have proved that, when T possesses no symmetries (i.e., has no non-trivial automorphisms), then k ∈ fs(T ) for all k ≥ xs(T ) + 1. This result is based on the explicit design of a distributed protocol enabling perpetual graph searching by k robots, for any k ≥ xs(T ) + 1.

When T possesses symmetries, the computation of fs(T ) becomes more complex. In this case, we have shown that fs(T ) depends on the set S T of isomorphic nodes of T separated by a path of even length. In particular, k / ∈ fs(T ) for any k ∈ [|S T |, |S T | + xs(T 0 )], where T 0 results from T after removing all nodes appearing in S T . Our impossibility results are not based on the perpetual nature of graph searching in min-CORDA, but on the exclusivity property. For the simpler case of a line, we have fully characterized fs(L) for any n-node line L. All our positive results are constructive.

We note that exclusive graph searching behaves very different from the classical one. As a result and due to min-CORDA, the proofs and the algorithms we propose are very different and more involved than in the classical case. Our proofs introduce several techniques that may prove useful also in future studies of exclusive graph-searching and min-CORDA.

Additional support from the ANR project "SHAMAN". ⋆ ⋆ ⋆ Partially supported by FP7 STREP EULER (N.N.)


Parikshit Dutta
email: parikshit.dutta@inria.fr.

Abhishek Halder
email: ahalder@tamu.edu.

Raktim Bhattacharya
email: raktim@tamu.edu.

Uncertainty Quantification for Stochastic Nonlinear Systems using Perron-Frobenius Operator and Karhunen-Loève Expansion

In this paper, a methodology for propagation of uncertainty in stochastic nonlinear dynamical systems is investigated. The process noise is approximated using Karhunen-Loève (KL) expansion. Perron-Frobenius (PF) operator is used to predict the evolution of uncertainty. A multivariate Kolmogorov-Smirnov test is used to verify the proposed framework. The method is applied to predict uncertainty evolution in a Duffing oscillator and a Vanderpol's oscillator. It is observed that the solution of the approximated stochastic dynamics converges to the true solution in distribution. Finally, the proposed methodology is combined with Bayesian inference to estimate states of a nonlinear dynamical system, and its performance is compared with particle filter. The proposed estimator was found to be computationally superior than the particle filter.

I. INTRODUCTION

Consider the stochastic nonlinear system dx (t) = f (x (t) , δ) dt + g w (x (t) , δ) dw (t) , (1) y (t) = h (x (t) , δ) + g η (x (t) , δ) η (t) ,

where x (t) ∈ R n and y (t) ∈ R m are the state and measurement vectors at time t, δ ∈ R p is the parameter vector, and w (t) ∈ R q , η (t) ∈ R l are mutually independent Wiener processes denoting process and measurement noise, respectively. The functions f (., .), h (., .) represent the dynamics and measurement model while the function-valued matrices g w (., .), g η (., .) denote respective noise coupling. If f (., .) and h (., .) are linear with constant noise coupling, the uncertainty quantification reduces to updating the mean and covariance propagation equations, and the associated estimation problem is solved by the well-known Kalman filter [START_REF] Kalman | A New Approach to Linear Filtering and Prediction Problems[END_REF]. For general nonlinear systems, if the uncertainty in the initial conditions (x 0 ) and parameters (δ) are prescribed through the joint probability density function (PDF) ρ 0 (z), z := [x δ]

T ∈ R n+p , then the propagation of uncertainty through (1) and (2) requires solving the Fokker-Planck equation (or Kolmogorov forward equation)

∂ρ ∂t = - n i=1 ∂ ∂x i (ρf i ) + n i=1 n j=1 ∂ 2 ∂x i ∂x j g w Qg T w ij ρ (3)
to compute ρ (z, t) subject to ρ (z, 0) = ρ 0 (z), where E [dw i ] = 0 and E [dw i dw j ] = Q ij 0, ∀ i, j = 1, 2, . . . , q.

Solving the parabolic partial differential equation (PDE) [START_REF] Bellman | Dynamic Programming[END_REF], second order in space and first order in time, is computationally hard [START_REF] Kumar | A Semianalytic Meshless Approach to the Transient Fokker-Planck Equation[END_REF] and the approximation algorithms often suffer from the 'curse of dimensionality' [START_REF] Bellman | Dynamic Programming[END_REF]. However, in the absence of process noise, (3) reduces to the Liouville equation associated with the Perron-Frobenius operator [START_REF] Lasota | Chaos, Fractals and Noise: Stochastic Aspects of Dynamics[END_REF] that describes the spatio-temporal drift of ρ (z, t) in the presence of parametric and initial condition uncertainty. Being a first order quasi-linear PDE, Liouville equation can be solved exactly [START_REF] Halder | Dispersion Analysis in Hypersonic Flight During Planetary Entry Using Stochastic Liouville Equation[END_REF] using the method of characteristics (MOC) [START_REF] Evans | Partial Differential Equations[END_REF] and unlike brute-force Monte-Carlo, such a computation does not compromise with runtime complexity [START_REF] Halder | Performance Bounds for Dispersion Analysis: A Comparison Between Monte Carlo and Perron-Frobenius Operator Approach[END_REF], [START_REF] Halder | Beyond Monte Carlo: A Computational Framework for Uncertainty Propagation in Planetary Entry, Descent and Landing[END_REF]. In nonlinear estimation setting, the prior PDF obtained through this propagation, can be used to compute the posterior PDF via Bayesian update. The resulting nonlinear filter has been shown [START_REF] Dutta | Hypersonic State Estimation using Frobenius-Perron Operator[END_REF] to outperform particle filter and bootstrap filter.

Instead of nonparametric propagation of the PDF, one may use a polynomial chaos [START_REF] Soize | Physical Systems with Random Uncertainties: Chaos Representations with Arbitrary Probability Measure[END_REF] expansion to quantify uncertainty, where a set of deterministic ordinary differential equations (ODEs) is derived using either Galerkin projection [START_REF] Spanos | Stochastic Finite Element Expansion for Random Media[END_REF] or stochastic collocation [START_REF] Xiu | High-Order Collocation Methods for Differential Equations with Random Inputs[END_REF]. There are three bottlenecks with this parametric approach. First, it's difficult to derive the deterministic ODEs for general nonlinearities [START_REF] Debusschere | Numerical Challenges in the Use of Polynomial Chaos Representations for Stochastic Processes[END_REF]. Second, the finite dimensional approximation of the probability space degrades the accuracy of long-term statistics. Third, the deterministic set of ODEs requires solution in a higher dimensional state space than the original, making the computation prohibitively expensive for large scale systems. Other parametric methods like moment closure lack accuracy due to approximation of the PDF using first few moments.

The objective of this paper is to introduce a mixed parametric-nonparametric approach for propagating uncertainty in the presence of process noise. Starting from the Langevin equation [START_REF] Langevin | On the Theory of Brownian Motion[END_REF], by N t -term Karhunen-Loève expansion [START_REF] Karhunen | Über lineare Methoden in der Wahrscheinlichkeitsrechnung[END_REF] of the process noise, we derive a deterministic ODE corresponding to the Itô stochastic differential equation (SDE) [START_REF] Kalman | A New Approach to Linear Filtering and Prediction Problems[END_REF]. Then the Perron-Frobenius operator based formulation is used to exactly solve the approximated SDE in distributional sense. A detailed verification procedure is illustrated to assess the quality of the N t -term approximation. There are two main contributions:

1) We show that for the nonlinear stochastic system (1) and ( 2), instead of solving (3), ρ (z, t) can be computed along the trajectories of an approximately equivalent deterministic dynamics over the extended state space. This enables us to harness the computational benefits of the Perron-Frobenius operator based formulation even in the presence of process noise. 2) We present a method for probabilistic verification of the approximated solution based on multivariate hypothesis testing. Consequently, for a fixed number of terms in the Karhunen-Loève expansion, we can compute the significance level of the finite-sample computation. Conversely, if a desired level of verification is specified, the number of terms for approximating the process noise can be computed. This paper is structured as follows. In section II, we show that the nonparametric propagation of the joint PDF ρ (z, t) can be achieved by solving the Liouville equation associated with the Perron-Frobenius operator, subject to the deterministic dynamics obtained from the Langevin equation by approximating the process noise through Karhunen-Loève expansion. Some examples are given in section III to illustrate the proposed formulation. Section IV describes the verification framework followed by applications to nonlinear filtering in section V. Section VI concludes the paper.

II. METHODOLOGY

We consider the Langevin equation

d dt x (t) = f (x (t) , δ) + ξ (t) , (4)
associated with the Itô SDE (1) for purely additive Gaussian white noise ξ (t) with autocorrelation σ 2 I. In the sequel, we first describe the method for uncertainty propagation via Perron-Frobenius operator followed by Karhunen-Loève expansion of the process noise. Next, we combine these two to compute the evolution of ρ (z, t) corresponding to (4).

A. Perron-Frobenius Operator

Definition 1: (Continuous-time Perron-Frobenius operator) Let X be endowed with a measure space (X , B, µ) and that all transformations f t of a semi-dynamical system [START_REF] Lasota | Chaos, Fractals and Noise: Stochastic Aspects of Dynamics[END_REF] {f t } t 0 are nonsingular (i.e. µ f -1 t (B) = 0, ∀B ∈ B such that µ (B) = 0), then for each t 0, the Perron-Frobenius operator P t : L 1 (X ) → L 1 (X ), corresponding to the dynamics f t , is uniquely defined as

B P t ρ (x) µ (dx) = f -1 t (B) ρ (x) µ (dx) , for B ∈ B. ( 5
)
The following properties of the Perron-Frobenius operator are important in our context, ∀ρ,

ρ 1 , ρ 2 ∈ L 1 , α, β ∈ R. 1) (Linearity) P t (αρ 1 + βρ 2 ) = αP t ρ 1 + βP t ρ 2 ,
2) (Non-negativity) P t ρ 0, provided ρ 0;

3) (Conservation)

X P t ρ (x) µ (dx) = X ρ (x) µ (dx) ,
The above properties guarantee that the operator

P t : L 1 (X ) → L 1 (X )
, is a Markov operator. In particular, property 3 assures that during Markov evolution, PDF respects the normality condition

X ρ (x) µ (dx) = 1
, where µ (.) is the Lebesgue measure.

One can show [START_REF] Lasota | Chaos, Fractals and Noise: Stochastic Aspects of Dynamics[END_REF] that the spatio-temporal evolution of the joint PDF via Perron-Frobenius operator [START_REF] Halder | Dispersion Analysis in Hypersonic Flight During Planetary Entry Using Stochastic Liouville Equation[END_REF], subject to the dynamics (4) with ξ (t) = 0, satisfies the Liouville PDE

∂ρ(x, t) ∂t = - n i=1 ∂ ∂x i (ρ(x, t)f i ) , (6)
which is simply (3) except the diffusion term. Being a first order PDE, ( 6) can be exactly solved using MOC along the characteristic curves. It can be shown [START_REF] Halder | Dispersion Analysis in Hypersonic Flight During Planetary Entry Using Stochastic Liouville Equation[END_REF] that the characteristic curves for ( 6) are nothing but the trajectories of the extended dynamics ż = [f (z) 0 p×1 ] T . The computational performance of this formalism is known [START_REF] Halder | Performance Bounds for Dispersion Analysis: A Comparison Between Monte Carlo and Perron-Frobenius Operator Approach[END_REF] to be superior than Monte-Carlo. However, in its original form, only initial condition and parametric uncertainties can be accounted for.

B. Karhunen-Loève Expansion of Process Noise

Karhunen-Loève (KL) expansion was derived independently by many researchers [START_REF] Karhunen | Über lineare Methoden in der Wahrscheinlichkeitsrechnung[END_REF], [START_REF] Loève | Fonctions Aléatoires du Second Ordre[END_REF] to represent a stochastic process Y (ω, t) as a random linear combination of a set of orthonormal deterministic

L 2 functions {e i (t)} ∞ i=1 , i.e. Y (ω, t) = ∞ i=1 Z i (ω) e i (t) . (7)
The idea is similar to the Fourier series expansion, where a deterministic linear combination of orthonormal L 2 functions is used. It's worth mentioning that in [START_REF] Halder | Performance Bounds for Dispersion Analysis: A Comparison Between Monte Carlo and Perron-Frobenius Operator Approach[END_REF], one can replace t by x ∈ R n , to have KL expansion of a random field [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF]. Further, if we write

Z i (ω) = √ λ i ζ i (ω), where λ i ∈ R + and {ζ i (ω)} ∞
i=1 is a sequence of random variables to be determined, then {λ i } ∞ i=1 and {e i (t)} ∞ i=1 can be interpreted as the eigenvalues and eigenfunctions of the covariance function [START_REF] Cressie | Statistics for Spatial Data[END_REF]

C (t 1 , t 2 ) := cov (Y (ω, t 1 ) -E [Y (ω, t 1 )] , Y (ω, t 2 ) -E [Y (ω, t 2 )]),
that admits a spectral decomposition [START_REF] Grigoriu | Stochastic Claculus: Applications in Science and Engineering[END_REF] of the form

C (t 1 , t 2 ) = ∞ i=1
λ i e i (t 1 ) e i (t 2 ). Since the covariance function is bounded, symmetric and positive-definite, the eigenvalue problem can be cast as a homogeneous Fredholm integral equation of second kind, given by,

Dt C (t 1 , t 2 ) e i (t 1 ) dt 1 = λ i e i (t 2 ).
Given the covariance function of a stochastic process, the eigenvalue-eigenfunction set can be found by solving Fredholm equation, and the resulting expansion [START_REF] Halder | Performance Bounds for Dispersion Analysis: A Comparison Between Monte Carlo and Perron-Frobenius Operator Approach[END_REF] converges to Y (ω, t) in mean-square sense [START_REF] Leon-Garcia | Probability, Statistics and Random Processes for Electrical Engineering[END_REF]. In particular, the following result is useful for our purpose.

Theorem 1: (KL expansion of Wiener process) (p. 548, [START_REF] Leon-Garcia | Probability, Statistics and Random Processes for Electrical Engineering[END_REF]) For Wiener process, the eigenvalues and eigenfunctions of the covariance function

C (t 1 , t 2 ) = (t 1 ∧ t 2 ), t 1 , t 2 ∈ [0, T ] × [0, T ], is given by λ i = 4T π 2 (2i -1) 2 , e i (t) = √ 2 sin i - 1 2 πt T , (8)
for i = 1, 2, . . ., and hence the KL expansion [START_REF] Halder | Performance Bounds for Dispersion Analysis: A Comparison Between Monte Carlo and Perron-Frobenius Operator Approach[END_REF] for w (t), is of the form

w (ω, t) m.s. = √ 2 ∞ i=1 ζ i (ω) sin i - 1 2 πt T i - 1 2 π T , (9)
where ζ i (ω) are i.i.d. samples drawn from N 0, σ 2 . Corollary 2: (KL expansion of Gaussian white noise) Since dw (t) = ξ (t) dt, the KL expansion for Gaussian white noise ξ (t) can be obtained by taking the derivative of ( 9) with respect to t, i.e.

ξ (t) m.s. = √ 2 ∞ i=1 ζ i (ω) cos i - 1 2 πt T , (10)
where the i.i.d. random variables

ζ i (ω) ∼ N 0, σ 2 .
By substituting [START_REF] Soize | Physical Systems with Random Uncertainties: Chaos Representations with Arbitrary Probability Measure[END_REF] in ( 4), we obtain a deterministic approximation of the stochastic dynamics (1). More generally, expansions of this nature can be derived in closed form, when the stochastic forcing in (1) can be written as a functional of Wiener process (e.g. Brownian bridge [START_REF] Karatzas | Brownian Motion and Stochastic Calculus[END_REF]).

C. KLPF Formulation

Let us consider the dynamical system given in (4) where x ∈ R n are the states, and z ∈ D z is the vector of state and parametric uncertainty, having initial PDF ρ 0 (z). Without loss of generality, it is assumed, that the system has only initial state uncertainty, hence z(t = 0) := x(t = 0). The stochastic forcing term ξ(t) in ( 4) is approximated using KL expansion having finite number of terms (N t ). In the present case, where ξ(t) is a Gaussian white noise, the approximation

reduces to, ξ (t) = √ 2 Nt i=1 ζ i (ω) cos i - 1 2 πt T , where t ∈ [0, T ].
After substituting this expression in (4), we get the following differential equation,

ẋ (t) = f (x (t)) + √ 2 Nt i=1 ζ i (ω) cos i - 1 2 πt T , (11)
Using the methodology given in [START_REF] Halder | Dispersion Analysis in Hypersonic Flight During Planetary Entry Using Stochastic Liouville Equation[END_REF], MOC is applied to [START_REF] Evans | Partial Differential Equations[END_REF], and an augmented dynamical system is formed, with states [x 1 (t) , . . . , x n (t) , ρ (x (t))] T , where ρ (x (t)) is the probability density of the states at time t, and x (t) = [x 1 (t) , . . . , x n (t)] T . The augmented system is given by,

ẋ (t) = f (x (t)) + √ 2 Nt i=1 ζ i (ω) cos πt(2i -1) 2T , ρ(x (t)) = -div f (x (t)) ρ(x (t)), (12)
where

div f (x (t)) = n i=1 ∂f (x (t)) ∂x i (t)
. Equation ( 12) can be solved to get the value of ρ (x (t)) along the characteristic curves of [START_REF] Evans | Partial Differential Equations[END_REF]. Detailed discussion of the solution methodology has been omitted here, and can be found in [START_REF] Halder | Beyond Monte Carlo: A Computational Framework for Uncertainty Propagation in Planetary Entry, Descent and Landing[END_REF], [START_REF] Dutta | Hypersonic State Estimation using Frobenius-Perron Operator[END_REF].

III. ILLUSTRATIVE EXAMPLES

The proposed methodology is applied to a Vanderpol's oscillator and a Duffing oscillator. The augmented dynamical system for the Vanderpol's oscillator is given by,

ẋ1 (t) = x 2 (t) (13a) ẋ2 (t) = 1 -x 2 1 (t) x 2 (t) -x 1 (t) + ξ (t) (13b) ρ (x (t)) = -1 -x 2 1 (t) ρ (x (t)) (13c)
and for the Duffing oscillator, is given by,

ẋ1 (t) = x 2 (t) (14a) ẋ2 (t) = 10x 1 (t) -30x 3 1 (t) -10x 2 (t) + ξ (t) (14b) ρ (x (t)) = 10ρ (x (t)) (14c)
with ξ (t) having autocorrelation 2πI. The initial state uncertainty, has a PDF ρ 0 (z) ∼ N [0, 0] T , diag (1, 1) for both the systems. Next, the initial PDF, ρ 0 (z) is sampled, with sample size of N = 5000. For the Vanderpol's oscillator, final time T is 1s, for the Duffing oscillator T = 3s. Number of terms in the KL expansion is fixed to N t = 7. Figure [START_REF] Kalman | A New Approach to Linear Filtering and Prediction Problems[END_REF] shows the evolution of probability densities with time for the two oscillators. The density value ρ (x (t)) is color coded with red representing high density value and blue representing low. It is observed that, for the Vanderpol's oscillator the probability mass accumulates along the limit cycle and for the Duffing oscillator, we get a bimodal PDF at final time. This is in agreement with the physical intuition of behavior of these systems.

IV. VERIFICATION OF SOLUTION

In the proposed methodology, the process noise in (4) has been approximated using finite number of terms in KL expansion. Due to this fact, the solution obtained after propagation is not error-free. It is well known that, KL expansion converges in mean square sense, to the underlying stochastic process as N t → ∞. But, same argument cannot be extended regarding the convergence of states x (t) of the dynamical system in [START_REF] Spanos | Stochastic Finite Element Expansion for Random Media[END_REF]. To obtain a meaningful solution, there should be some notion of convergence of the states in [START_REF] Spanos | Stochastic Finite Element Expansion for Random Media[END_REF] to states in (4). Hence, verification of the obtained solution is important in this scenario. In the present work, we have verified that, given the finite term KL expansion of ξ (t), the solution of the approximated dynamics in [START_REF] Spanos | Stochastic Finite Element Expansion for Random Media[END_REF] converge to the true solution in distribution. We have used the Kolmogorov-Smirnov (KS) test to verify the solution of approximated dynamics, as detailed below.

The Kolmogorov-Smirnov test [START_REF] Justel | A Multivariate Kolmogorov-Smirnov Test of Goodness of Fit[END_REF] is a statistical test used to compare a sample with a reference probability distribution. It quantifies a distance D N , between reference cumulative distribution function (CDF), F (x), and empirical CDF of the sample being tested, F N (x), which is given by D N = sup x |F N (x) -F (x)|. Here, N refers to sample size of the given sample. The null hypothesis is that the sample comes from the reference distribution. Given a significance level α, the null hypothesis is accepted if Here, K is a random variable that follows Kolmogorov distribution, with CDF

√ N D N ≤ K α , where, Pr (K ≤ K α ) = 1 -α. (15)
Pr (K ≤ x) = √ 2π x ∞ i=1 e -(2i-1) 2 π 2 8x 2 . ( 16
)
If the null hypothesis is accepted, then the KS test is passed and it can be concluded that the given sample comes from the reference PDF.

Due to the fact that N < ∞, D N is a random variable, as each sample will give us different

D N values. The empirical distribution of √ N D N is given by, F Ns (x) = 1/N s ( Ns i=1 I x< √ N D i N
), where N s is the number of samples and D i N is the D N value for i th sample. The Dvoretzky-Kiefer-Wolfowitz inequality [START_REF] Dvoretzky | Asymptotic Minimax Character of the Sample Distribution Function and of the Classical Multinomial Estimator[END_REF] characterizes the rate-ofconvergence of F Ns (x) as [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF] where F K (x) is the Kolmogorov CDF in [START_REF] Loève | Fonctions Aléatoires du Second Ordre[END_REF]. Hence as N s → ∞ the test statistic √ N D N exponentially converges in distribution to a Kolmogorov random variable with rate 2e -2Nsε 2 .

Pr sup x∈R |F Ns (x) -F K (x)| > ε ≤ 2e -2Nsε 2 , ∀ε > 0,
In the present case, the analytical representation of PDF at time t = T is not known. Hence to verify our solution, we back-propagate the sample obtained after propagation of (11) at final time T to time t = 0 using the original dynamics in (4); and check if the back-propagated sample belongs to the initial PDF ρ 0 (z).

Let the elements of the sample; sampled from the initial PDF be x 0,i , i = 1, . . . , N , and elements after propagation of (11) be x T,i . Let the back-propagated sample obtained by propagating (4) from [T, 0] be x0,i . Here, we claim that, x0,i is a sample from the initial PDF. Hence, using the KS test, the sample x0,i , is compared with the initial PDF, for statistical significance. Here, we have used a multivariable KS test, which is similar to the one presented in [START_REF] Justel | A Multivariate Kolmogorov-Smirnov Test of Goodness of Fit[END_REF]. The verification methodology is described in algorithm 1.

We apply the proposed verification methodology to the Vanderpol's oscillator given in [START_REF] Debusschere | Numerical Challenges in the Use of Polynomial Chaos Representations for Stochastic Processes[END_REF]. The initial state uncertainty has the PDF ρ 0 (z) ∼ N ([0, 0], diag (1, 1)). We obtain N s = 100 samples, each of sample size N = 500 from ρ 0 (z). The number of terms in the KL expansion is fixed to N t = 21. At first, we pick a sample from the 100 available Reject H0 for the j th sample 16:

end if 17: end for ⊲ Repeat the same test for the next sample samples and apply the verification methodology. Figure [START_REF] Kumar | A Semianalytic Meshless Approach to the Transient Fokker-Planck Equation[END_REF] shows plot for the location of elements initially (x 0 ) and after back-propagation (x 0 ), for the particular sample. It is observed that the back-propagated sample, is clustered around the origin and sparseness increases as we move away. This is in agreement with the physical intuition of random samples drawn from ρ 0 (z) which is a standard normal distribution.

Using the methodology in [START_REF] Justel | A Multivariate Kolmogorov-Smirnov Test of Goodness of Fit[END_REF] we obtain uniformly distributed sample ŷ0 = F (x 0 ) using inverse transform sampling theory [START_REF] Steinbrecher | Quantile mechanics[END_REF]. Figure [START_REF] Kumar | A Semianalytic Meshless Approach to the Transient Fokker-Planck Equation[END_REF] shows the plot of empirical CDF of ŷj 0 and the uniform CDF for the randomly selected sample. It can be seen that, visually the CDFs are close to each other. Hence, the value of D N for the particular sample is expected to be less. For the KS test, the value of α is fixed to 0.05. Hence, using ( 15) and ( 16), we get Kα √ N = 0.0607. For the given sample, the value of D N was 0.0455. Hence, the null hypothesis, that the sample is drawn from a standard normal distribution, is accepted with a significance level α = 0.05.

Figure [START_REF] Bellman | Dynamic Programming[END_REF] show plot for D N values for all the 100 samples drawn. The samples above red line fail the test if α = 0.05 and the samples above the black line fail, if α = 0.01. It is observed that only 28 samples fail the test when α = 0.05 and 14 samples fail when α = 0.01. It can be seen that, as the significance level α, is decreased, the chances that the null hypothesis is accepted increases. As mentioned before in [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF], the value of D N obtained after KS test is a random variable, whose CDF converges to Kolmogorov distribution exponentially. we do not observe a monotonic behavior. But, there are very few outliers, as the probability in the left hand side of ( 17) can be decreased arbitrarily by varying ǫ.

V. APPLICATION TO NONLINEAR FILTERING

We apply the proposed KLPF methodology to estimate states of a nonlinear system. The nonlinear estimation algorithm used is same as the one described in [START_REF] Dutta | Hypersonic State Estimation using Frobenius-Perron Operator[END_REF]. Here, we draw random particles (elements) from the domain D z of z(t = 0). Propagation of uncertainty is done using [START_REF] Xiu | High-Order Collocation Methods for Differential Equations with Random Inputs[END_REF]. By using this framework for nonlinear estimation purposes, we get the exact value of prior PDF after propagation. This gives us significant advantage over particle filters where, generally the prior PDF is approximated using histograms. Bayesian inference is used to update the posterior PDF from prior for each particle.

In the present work, the system used is the Vanderpol's oscillator described by [START_REF] Debusschere | Numerical Challenges in the Use of Polynomial Chaos Representations for Stochastic Processes[END_REF]. We use a nonlinear measurement model given by, y (t) = x 2 1 (t) + x 2 2 (t) + η (t), where η (t) is zero mean Gaussian measurement noise with autocorrelation R = 6 × 10 -2 . The process noise, ξ (t) is also considered to be Gaussian with zero mean and an autocorrelation Q = 6 × 10 -1 .

We consider the initial state uncertainty to be normally distributed with E[z] = [0, 0] T and E[z 2 ] = diag [START_REF] Kalman | A New Approach to Linear Filtering and Prediction Problems[END_REF][START_REF] Kalman | A New Approach to Linear Filtering and Prediction Problems[END_REF]. It is assumed that the initial states of the actual system is [0.25, 0.25], with initial error in estimation being [0.25, 0.25]. The measurement update interval was fixed to 0.1s and final time T was assumed to be T = 1s.

The performance of the KLPF-based filter is compared with the generic particle filter given in [START_REF] Arulampalam | A Tutorial on Particle Filters for Online Nonlinear/Non-Gaussian Bayesian Tracking[END_REF]. The sample size for each estimator is fixed to 5000 elements. The simulations

Fig. 1 .

 1 Fig. 1. Uncertainty Propagation for Vanderpol's Oscillator (13) in top row, and Duffing Oscillator (14) in bottom row.

Algorithm 1 7 :

 17 Verification of KL expansionRequire: Domain Dz pdf ρ0(z) of initial parameters. Significance level α, sample size N , Number of samples Ns, Number of terms in KL expansion Nt. 1: Calculate Kα using (15) 2: for j = 1 to Ns do ⊲ Sample number Counter 3: Draw N elements x j 0,i at random from Dz. ⊲ Use Markov Chain Monte Carlo 4: Create null hypothesis H0 : ρ = ρ0 for the current sample 5: for i = 1 to N do ⊲ Sample element counter 6: Propagate x j 0,i , from t = [0, T ], using (11) to get x j T,i ⊲ Use of KL approximated dynamics with number of terms Nt Backward propagate x j T,i , from t = [T, 0] using (4) to get xj 0,i . ⊲ Use of original Stochastic Dynamics 8:Get ŷj 0,i ← ρ0(x j 0,i ) ⊲ From inverse transform sampling theory, ŷj 0 is uniformly distributed in [0, 1]

Fig. 2 .

 2 Fig. 2. (Top) Scatter plot of initial and the back-propagated elements for a given sample. Blue circles are the original sample and the red ones are the back-propagated sample. (Bottom) Plot of Empirical CDF of ŷj 0 (left) and uniform CDF (right) for a given sample

Fig. 3 .

 3 Fig. 3. D N value for all the 100 samples. The samples above the red and black lines fail when α = 0.05 and 0.01 respectively.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Empirical CDF of √ N D N (red) and Kolmogorov CDF (blue)

This research was partially supported by NSF award # 1016299 with D. Helen Gill as the Program Manager.

were performed on a Linux Machine with Intel R Pentium D processor. Figure [START_REF] Evans | Partial Differential Equations[END_REF] and Figure [START_REF] Halder | Performance Bounds for Dispersion Analysis: A Comparison Between Monte Carlo and Perron-Frobenius Operator Approach[END_REF] show the ±3σ plots for the estimators. No major difference can be observed in the performance of the estimators as the errors in estimation are within the ±3σ limits and converge. But if we compare the computational time, KLPF-based estimator takes 58.46s per filtering step, whereas the time taken by generic particle filter is 284.19s. Thus it can be concluded that, given same sample size and the same processing environment, KLPFbased estimator achieves similar performance as the generic particle filter, and is almost five times faster than the particle filter based estimator.

VI. CONCLUSIONS AND FUTURE WORK

In this work, an uncertainty propagation methodology is proposed for stochastic nonlinear dynamical systems, based on Karhunen-Loève (KL) expansion and Perron-Frobenius operator. The methodology uses KL expansion to approximate the process noise and Perron-Frobenius operator to predict the evolution of uncertainty. The method is illustrated by propagating uncertainty in a Duffing oscillator and a Vanderpol's oscillator. It was observed that the probability distribution of states, obtained after propagation, is in compliance with the underlying theory for these nonlinear systems. Moreover, the solution obtained with the approximated dynamics using KL expansion, was verified with the true solution, when applied to a Vanderpol's oscillator. It was found that the states of the approximated system converges in distribution to the true states of the system. A verification algorithm for the same has been presented here.

In future, we will focus on obtaining convergence relation between approximated dynamics and the true solution as a function of the number of terms in KL expansion.