A Sieve Method for Consensus-typeNetwork TomographyThe research of the authors was supported by AFOSR grant FA9550-09-1-0091.

# A Sieve Method for Consensus-type Network Tomography††thanks: The research of the authors was supported by AFOSR grant FA9550-09-1-0091.

Marzieh Nabi-Abdolyousefi and Mehran Mesbahi The authors are with the Department of Aeronautics and Astronautics, University of Washington, Seattle, WA 98195-2400 USA (emails: mnabi+mesbahi@uw.edu)
###### Abstract

In this note, we examine the problem of identifying the interaction geometry among a known number of agents, adopting a consensus-type algorithm for their coordination. The proposed identification process is facilitated by introducing “ports” for stimulating a subset of network vertices via an appropriately defined interface and observing the network’s response at another set of vertices. It is first noted that under the assumption of controllability and observability of corresponding steered-and-observed network, the proposed procedure identifies a number of important features of the network using the spectrum of the graph Laplacian. We then proceed to use degree-based graph reconstruction methods to propose a sieve method for further characterization of the underlying network. An example demonstrates the application of the proposed method.

Keywords: Inverse problems, coordination algorithms, system identification, graph reconstruction.

## I Introduction

Physical sciences are often concerned with inferring models and physical parameters from data. Given a model for a physical phenomena, computing the data values is often referred to as the forward problem. On the other hand, in inverse problems, the objective is the construction, validation, invalidation, or reconstruction of the model from a set of measurements associated with the system. Inverse problems arise in fields such as astronomy, geophysics, medical imaging, remote sensing, ocean acoustic tomography, and non-destructive testing [1, 2, 3]. Closer to the present work are the inverse problems associated with electrical networks [4], and the celebrated “Can one hear the shape of a drum?" which aims to characterize a manifold via its spectra [5], or more recently, “Can one hear the shape of a graph?" [6]. In fact, in this paper, we address the inverse problem related to consensus-type coordination algorithms. Consensus-type algorithms have recently been employed for analysis and synthesis of a host of distributed protocols and control strategies in multi-agent systems, including, flocking, formation control, rendezvous, and distributed estimation [7].

One of the key aspects of this class of protocols is the strong dependency between the interaction and information-exchange geometry among the multiple agents, on one hand, and the dynamic properties that these systems exhibit, on the other. Motivated by this dependency, in our work, we consider the scenario where the interaction network is inside a “black box,” and that only certain “boundary” nodes in the network can be influenced and subsequently observed. The “input” boundary nodes are then used to stimulate the network, whose response is subsequently observed at the “output” boundary nodes. Using this setup, in our complementary work [8], we have presented a node knockout procedure that aims to find the generating function of the graph Laplacian from the observed input-output data. Our focus in the present work, in the meantime, is to reduce the search space for the identification of the network topology by blending ideas from system identification, integer partitioning, and degree-based graph reconstruction. The implicit contribution of our analysis is its ramifications for exact identification from boundary nodes for networks that have an embedded consensus-type algorithms for their operation, including formation flying, distributed estimation, and mobile robotics.

Our notation and terminology are standard.111The main focus of this work is on undirected graphs. However, the extension of some results to the directed case will be examined in subsequent works. We denote by the undirected simple graph with vertex set and edge set , comprised of two-element subsets of ; we use “nodes” or “agents” interchangeably with “vertices.” Two vertices are called adjacent if . For vertex , denotes the number of its adjacent vertices or neighbors. The Laplacian matrix for the graph is denoted by . Laplacian matrices are positive semi-definite whose spectrum will be ordered as . We use to denote the characteristic polynomial of the graph Laplacian. The cardinality of the set will be denoted by ; on the other hand denotes a function of that is bounded by some constant multiple of for large values of .

## Ii Network Identification

Consider the consensus protocol adopted by -nodes, where is the state of the -th node, e.g., its position, speed, heading, voltage, etc., evolves according to the sum of the differences between the -th node’s state and its neighbors. Next, let a group of agents with cardinality , “excite” the underlying coordination protocol by injecting signals to the network, with another set of agents , of cardinality , measuring the corresponding network response. Hence, the original consensus protocol from node ’s perspective assumes the form

 ˙xi(t)=∑{i,j}∈E(xj(t)−xi(t))+Biui(t), (1)

where if , and zero otherwise. Without loss of generality, we can always assume that and modify the control signal as if necessary. Adding the observation ports to this “steered” consensus, and having when , we arrive at the compact form of an input-output linear time-invariant system,

 ˙x(t)=A(G)x(t)+Bu(t),y(t)=Cx(t), (2)

where , , and .

Even though in general sets and can be distinct and contain more than one element, for the convenience of our presentation, we will assume that they are identical- and at times, assume that the resulting input-output system is in fact SISO. The extension of the presented results to the case when and are distinct will be discussed after introducing the basic setup and approach.

We now pose the inverse problem of graph-based coordination algorithms, namely, the feasibility of identifying the spectral and structural properties of the underlying network via the data facilitated by the input-output ports and . In order to implement this program, however, we need to assume that: (1) the identification procedure has knowledge of the number of agents in the network, and (2) the input/output sets and have been chosen such that the system described in (2) is controllable and observable. Although the first assumption is reasonable in general, the second one requires more justification which we now provide. In the trivial case when and is equal to the identity matrix, the input-output consensus (2) is controllable and by duality, observable. However, more generally, the controllability/observability of the network from a subset of its boundary nodes, is less trivial, and more to the point, not guaranteed for general graphs [7]. In the meantime, since we will need controllability and observability of the network for its identifiability, we will rely on a topical conjecture in the algebraic graph theory community to the effect that for large values of , the ratio of graphs with nodes that are not controllable from any single node to the total number of graphs on nodes approaches zero as  [9]; this phenomena is depicted in Fig. 1. In the present paper, we take the controllability and the observability of the underlying graph from the input and output nodes as our working assumption. In the meantime it is always convenient to know when the network is uncontrollable from a given node.

###### Lemma II.1.

Let as the input-output realization of (2). The uncontrollable/unobservable eigenvalues of (2) will not appear in the corresponding entry of . Specifically, will be order polynomial for the SISO case with agents and uncontrollable/unobservable eigenvalues.

###### Proof.

Since the underlying graph is undirected, matrix is symmetric and there exists a unitary matrix and a real nonnegative diagonal matrix such that . The columns of are an orthonormal set of eigenvectors for . The corresponding diagonal entries of are the eigenvalues of  [10]. Therefore,

 G(s)=C(sI−A(G))−1B=C(sI−UΛUT)−1B=CU(sI−Λ)−1UTB (3)

From PBH test, if the system (2) is not controllable, there is an eigenvector that is orthogonal to . Therefore for an arbitrary uncontrollable eigenvalue , the -th row of is orthogonal to and will not appear in . An analogous argument works for the unobservable case. ∎

### Ii-a System Identification

We now consider various standard system identification procedures in the context of identifying the spectra of the underlying graph Laplacian, and subsequently, gaining insights into the interconnection structure that underscores the agents’ coordinated behavior.

System identification methods are implemented via sampling of the system (2) at discrete time instances,222The system identification methods work based on data sampling from the system. Since we aimed to identify the interaction geometry of the network, we originally considered a continuous system. Therefore, we need to discretize the system (2). , with , assuming the form

where , , , , , and .333The notation for a square matrix refers to its matrix exponential. In fact, the system identification process leads to a realization of the model

where is the realization of in (4). The estimated system (5), on the other hand, is equivalent to the continuous-time system

 ˙˜x(t)=˜A˜x(t)+˜Bu(k),y(t)=˜C˜x(t), (6)

with and in this case, where denotes the matrix logarithm. Since the system (5) is a realization of the system (4), it follows that the estimated triplet is a realization of in (2). As a result, there exists a similarity transformation induced by the matrix , such that , , and . In fact, in the controllable/observable case, eigenvalues of are precisely matched with the eigenvalues of . Obtaining a zero as eigenvalue of , which is equivalent of obtaining as the eigenvalue of , is a sign of uncontrollable and/or unobservable mode in (2).444This follows from Lemma II.1 since will appear as zero in the corresponding entries of . For example in the identification procedure called Iterative Prediction-Error Minimization Method, the model (4) for every input and output can be represented as , where and . The unknown model parameters can then be estimated by comparing the actual output and predicted output using the mean-square minimization. In this case, the output predictor is constructed as . In yet another candidate system identification procedure, namely the Subspace Identification Method, the system (4) is approximated by another system in the form (5), using the state trajectory of the dynamic system that has been determined from input-output observations. The Hankel matrix, which can be constructed from the gathered input-output data, plays an important role in this method. By constructing the Hankel matrix, the discrete time system matrices , , and can then be determined. Subsequently, the continuous-time estimated matrices , , and can be identified; see [11] for an extensive treatment of system identification methods.

In summary, an identification procedure such as the above two methods, implemented on a controllable and observable steered-and-observed coordination protocol (2), leads to a system realization whose state matrix is similar to the underlying graph Laplacian and in particular sharing the same spectra and characteristic polynomial. However, a distinct and fundamental issue in our setup is that having found a matrix that is “similar” to the Laplacian of a network is far from having exact knowledge of the network structure itself [8]. This observation motivates the following question: to what extend does the knowledge of the spectra of the graph, combined with the knowledge of the input-output matrices, reduce the search space for the underlying interaction geometry? In this note, we explore this question using techniques based on integer partitioning and degree-based graph reconstruction.

## Iii Graph Characterization

We first review qualitative characterization of the underlying interconnection topology via its identified characteristic polynomial. We then explore the possibility of reducing the search space for the underlying network via the proposed sieve method.

### Iii-a Graph Characterization via the Characteristic Polynomial

Recall that via a system identification method, the characteristic equation of the system (2) can be found as

 ϕG(s) = det(sI−A(G))=sn+a1sn−1+...+an−1s+an. (7)

Although the spectra of the graph Laplacian in general is insufficient to form an explicit characterization of the underlying network, it leads to a number of useful structural information about its geometry; we list a few: (1) the value is the number of spanning trees in , (2) one has , where is the number of edges in the graph, (3) if , the underlying interconnection is a tree. For a tree, the coefficient is also called the graph Wiener index (the sum of all distances between distinct vertices of ) [12], (4) if the associated graph is a tree, is the number of -matching in the subdivision of (where each edge of is replaced by a path of length 2), (5) if the eigenvalues of are distinct, with where , define . It is then well-known that , where the graph is the complement of the graph . The Hoffman number of the graph, , can also be found as whose properties and applications have been studied in [13], (6) let be a tree with vertices. If has only one positive Laplacian eigenvalue with multiplicity one, then is the star  [13], (7) let be a connected graph with exactly three distinct Laplacian eigenvalues. Then the algebraic connectivity (the second smallest Laplacian eigenvalue) of is equal to one if and only if is a star of with , (8) if is a connected graph with integer Laplacian spectra, then , where is the diameter of the graph and .

Although the spectra of the Laplacian provides important insights, as noted above, into the structural properties of the network, we now proceed to explore the possibility of complete identification of the underlying network using its graph spectra complemented with a sieve method.

### Iii-B Graph Sieve

In this section, we provide an overview of the graph sieve procedure– that in conjunction with the identified Laplacian spectra– leads to a more confined search for the network in the black box. The essential ideas involve the judicious use of integer partitioning algorithms and degree-based graph reconstruction.

Recall that with the standing assumption of in (2), for the identified system matrices , after appropriate relabeling, the product leads to the first block partition of the matrix . Notice that if , we still obtain entries of the matrix which may not contain the diagonal entries. The product gives us the degree of the nodes that are in common in both input and output sets and information on the minimum degree of other nodes. Since the eigenvalues of the identified matrix are identical to those of , as the result of the identification process, we have access to the sum of the degrees of all nodes in the network, as well as the degrees of a subset of -boundary nodes. Let us define as the set of nodes in -boundary nodes which appears in both , and with . Moreover, let

 rd=∑v∈~R\bf deg vands=\bf trace(˜A)−rd. (8)

If , then , and the set will be equal to -boundary nodes. We can then proceed to determine the degrees of remaining nodes, or equivalently, partition the positive integer (8) into integers, each assuming a value between 1 and , and a lower bound for the degrees of nodes [14]. The possible values for the partitioning comes from the fact that we expect the resulting graph be connected while respecting the bounds on the maximum allowable node degrees.

Partitioning integers without constraints on the resulting partition is often referred to as unrestricted partitions. Restricted partitions, on the other hand, are those with constraints on the largest value of the partition that is no greater than a value of , or no smaller than , or both. Algorithms that generate unrestricted partitions can often be used to generate the restricted ones by certain modifications. Several such algorithms, dealing with unrestricted and restricted integer partitioning, have been suggested in literature. In the context of the graph realization using the proposed system identification method, we proceed to use the algorithms in [15], in order to generate different possible sets of integers between 1 and , and nodes with specified lower bounds on their degree such that their sum is (8).

### Iii-C Integer Partitioning Algorithms and Complexity Analysis

Consider a degree sequence while , with specified lower bound on of them. Without loss of generality, assume that the first degrees are lower bounded as

 di≥Lifori=1,…,r−~r. (9)

We are interested to find all possible partitioning of into integers between and satisfying (9). We have the following observation; see [16].

###### Lemma III.1.

Let the number of partitioning of into integers between and be denoted by . Then

 Pn−~r(s)=Pn−~r−1(s−1)+(n−~r)Pn−~r(s−1). (10)

The following algorithm, proposed in [15] finds all partitioning of into integers between and satisfying (9). The partitioning of with components can be generated in increasing lexicographic order by starting with and continuing as follows. To obtain the next partition from the first one, scan the elements from right to left, stopping at the right most such that . Replace by for and then replace by . For example, if we have and the partition , we find that is greater by than the rightmost , and so the next partition is . When no element of the partition differs from the last by more than , we are done.

In the suggested algorithm the output size of each partitioning of into some arbitrary number of integers , , is . This means that the total output size is . The approximate size of the number is provided by the following asymptotic formula,

 P(s)∼14√3exp(π√2s3).\vspace−0.10in

In other words, grows faster than any polynomial, but slower than any exponential function . However, in our application, we are interested in a subset of which has a specified size and satisfy certain constraints. Specifically, the integer in (8) is approximately the number of edges in the graph. For simple graphs if , then the upper bound for the proposed partitioning is , and if , the upper bound for the proposed partitioning is .

Although the algorithm above leads to a possible degree sequence for the underlying graph— consistent with the identification procedure– we need an additional set of conditions for ensuring that the obtained sequence in fact corresponds to that of a graph.

###### Definition III.2.

A graphical sequence is a list of nonnegative numbers that is the degree sequence of some simple graph. A simple graph with degree sequence is said to realize .

Our next step is therefore to characterize the necessary and sufficient conditions for a set of integers to be graphical. For this, we resort to the following result.

###### Theorem III.3.

[17] For , an integer list of size is graphical if and only if is graphical, where is obtained from by deleting its largest element and subtracting from its -th next largest elements. The only -element graphical sequence is

###### Example III.4.

Consider a sequence . Since the number of odd degree nodes is odd, the sequence is not graphical. Let us also construct the sequences described above, as and , which again verifies that this sequence is not graphical.

Since the maximum number of steps to check whether a sequence is graphical or not is , the complexity of this algorithm is . Given the particular algorithmic means of generating a graphical sequence for the required integer partitions, as detailed above, we now consider the problem of constructing graphs based on a graphical sequence.

### Iii-D Degree Based Graph Construction Algorithms and Complexity Analysis

Before describing the algorithm, we need to provide a few definitions.555The necessary definitions and algorithms have been discussed in [18] and are briefly described here to complement the presentation. Let denote the adjacency set of node defined as

 A(i)={ak|ak∈V,ak≥i,for allk, 1≤k≤di}.

The reduced degree sequence is obtained after removing node with all its edges from . We now define the ordering between two adjacency sets of node , and , as if we have for all . In this case we also say that is "to the left" of . The next lemma introduces a sufficient condition for the sequence to be graphical.

###### Lemma III.5.

[18] Let be a non-increasing graphical sequence, and let , be two adjacency sets for some node , such that . If the degree sequence reduced by (that is ) is graphical, then the degree sequence reduced by (that is ) is also graphical.

The above lemma guarantees preservation of “graphicality” for all adjacency sets to the left of a graphical one. Now consider a graphical degree sequence on nodes obtained from previously discussed integer partitioning approach. From the identity , as we discussed in § II, entries of the system matrix are known;666In the case where , entries of are known. define these set of edges as being “pre-determined” in the graph which cannot be repeated again. Put these connections in the forbidden set . The Algorithm 2 describes how we can construct all possible graphs avoiding the edges in .

When constructing , checking graphicality is only needed for those adjacency sets which are incomparable by the ordering relationship to any of the current elements of ; for the remaining sets graphicality is guaranteed by Lemma III.5.

The total number of graphs that the above algorithm produces is . Of course, this procedure makes sense for degree sequences for which there is only a small number of labeled graphs realizing it. An upper bound on the worst case complexity of the algorithm for constructing a sample from a given degree sequence is , with being the number of edges in the graph. For simple connected graphs, the maximum possible number of edges is , and the minimum possible number is . If , then , and if , then , which is an upper bound, independent of the degree sequence [19]. An algorithm to realize graphical degree sequences of directed graphs has been studied in [20] which can be use to extend these results to simple directed graphs.

Our sieve method for confining our search for the underlying graph following the system identification of § II, thus involves:

1. perform an integer partition on the value (8), keeping the partitions that lead to a graphical sequence; let denote the set of all graphs that remain after this first stage of the sieve.

2. construct a candidate connected graph in that is consistent with the matrix and satisfy the given degree sequence,

3. compare the Laplacian eigenvalues of the constructed graphs with the roots of the characteristic polynomial (7) and discard the inconsistent graphs.

The candidate graphs for the underlying network topology are now among the ones that remain after this three step sieve. We note that the sieve method is guaranteed to reduce the search space for the graph structure by at least a factor of , a bound that is obtained from the bound on the number of permissible degree-based integer partitioning facilitated by the network system identification. An example for this procedure is given next.

## Iv An Example

Our goal in this example is to gather information on the graph shown in Figure 2(a) using the system identification procedure. Using nodes , , and as the input-output nodes in (2), we obtain . Since the polynomial has just one zero root, the underlying graph is connected. Moreover, the graph is not a tree due to the fact that . The graph has edges and spanning trees. We also obtain the estimation matrices , and to be

 ˜A = ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣−2.09139.5928−3.2434−1.9321−8.3119−3.6056−0.7817−2.6237−2.3217−3.3284−2.10370.87470.2080−0.1114−3.47150.84430.0218−0.2108−0.3304−0.12690.5219−3.15000.09990.3707−0.64062.3385−3.4079−4.7957−7.18710.6395−1.87049.1263−2.8730−1.3239−7.4520−3.4763⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦ ˜B = ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣−9.790313.7604−4.02111.66500.0982−1.67520.63480.1419−0.31470.4699−0.3308−0.41160.34302.7623−3.1201−9.301912.8359−3.5526⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦ ˜C = ⎡⎢⎣−28.15422.38260.35470.88800.643430.0441−28.23112.00050.07290.27810.963430.1259−28.24552.5397−0.5074−0.50560.515730.1419⎤⎥⎦.

And also

 ˜C˜A˜B=⎡⎢⎣\par−3.00001.0000−0.00001.0000−4.00001.00000.00001.0000−3.0000⎤⎥⎦.

Since the diagonal of the matrix is , , , and . In the meantime, the sum of degrees of the remaining nodes is . The possible integer partitions for the remaining three nodes such that the sum is and each degree is less than will be , , and . Therefore, the set of the possible degrees sequences is comprised of , , and . According to the Theorem III.3, three set of integer partitioning are realization of some graphs.

Next, we construct the graphs with the three candidate degree sequences by implementing the algorithm in [18] and constructing all connected graphs consistent with , and with degree sequences , , and . The first degree sequence, , and the degree sequence are realizations of two graphs while the third degree sequence is a realization of two graphs. All these graphs satisfy the constraint imposed by . By comparing the Laplacian spectra for the corresponding graphs with the roots of the identified characteristic polynomial, we can thereby identify the original graph. Three candidates for the constructed graph are depicted in Fig. 2(b)-(d).

## V Conclusion

In this paper, we introduced a network identification scheme which involves the excitation and observation of nodes running consensus-type coordination protocols. Starting with the number of vertices in the network as a known parameter, as well as the controllability and observability of the resulting steered-and-observed network, the proposed procedure strives to collect pertinent information on the topology of the underlying graph. In this direction, we examined the applications of spectral characterization of graphs as well as a sieve method that is based on integer partitioning algorithms and feasible graphical sequences.

## References

• [1] S. Babaeizadeh, D. Brooks, and D. Isaacson, “3-d electrical impedance tomography for piecewise constant domains with known internal boundaries,” IEEE Transactions on Biomedical Engineering, vol. 54, no. 1, pp. 2–10, 2007.
• [2] O. Marques, T. Drummond, and D. Vasco, “A computational strategy for the solution of large linear inverse problems in geophysics,” in Parallel and Distributed Processing Symposium, April 2003, pp. 22–26.
• [3] J. Mueller, S. Siltanen, and D. Isaacson, “A direct reconstruction algorithm for electrical impedance tomography,” IEEE Transactions on Medical Imaging, vol. 21, no. 6, pp. 555–559, 2002.
• [4] E. B. Curtis and J. A. Morrow, Inverse problems for electrical networks.   USA: World Scientific, 2000.
• [5] M. Kac, “Can one hear the shape of a drum?” The American Mathematical Monthly, vol. 73, no. 4, 1966.
• [6] B. Gutkin and U. Smilansky, “Can one hear the shape of a graph?” Journal of Physics A: Mathematical and General, vol. 34, pp. 6061–6068, 2001.
• [7] M. Mesbahi and M. Egerstedt, Graph Theoretic Methods for Multi-Agent Networks.   Princeton University Press, Forthcoming.
• [8] M. Nabi-Abdolyousefi and M. Mesbahi, “Can one hear the shape of coordination?” 49th IEEE Conference on Decision, Atlanta, December 2010.
• [9] C. D. Godsil, “Controllablity on networks,” http://quoll.uwaterloo.ca/mine/Talks/control.pdf.
• [10] R. A. Horn and C. R. Johnson, Matrix Analysis.   Cambridge University Press, 1985.
• [11] L. Ljung, System identification.   NJ: Prentice-Hall: Upper Saddle River, 1999.
• [12] B. Mohar, “The laplacian spectrum of graphs,” Graph Theory, Combinatorics, and Applications, vol. 2, pp. 871–898, 1991.
• [13] Y. Teranishi, “The hoffman number of a graph,” Discrete Mathematics, vol. 260, no. 1-3, pp. 255–265, 2003.
• [14] L. E. Dickson, History of the theory of the numbers, Volume II, Diophantine Analysis.   New York: Chelsea Publishing Compnay, 1971.
• [15] E. M. Reingold, J. Nievergelt, and N. Deo, Combinatorial algorithms: theory and practice, 1977.
• [16] M. Bona, A walk Through Combinatorics.   World Scientific Publishing Co., 2002.
• [17] D. B. West, Introduction to Graph Theory.   Prentice Hall, 2001.
• [18] H. Kim, Z. Toroczkai, P. L. Erdos, I. Miklos, and L. A. Szekely, “Degree-based graph construction,” Journal of Physics A: Mathematical and Theoretical.
• [19] C. I. D. Genio, H. Kim, Z. Toroczkai, and K. E. Bassler, “Efficient and exact sampling of simple graphs with given arbitrary degree sequence,” arXiv:1002.2975v1, 2010.
• [20] P. L. Erdos, I. Miklos, and Z. Toroczkai, “A simple havel-hakimi type algorithm to realize graphical degree sequences of directed graphs,” arXiv:0905.4913v2, 2010.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters