A Sieve Method for Consensustype
Network Tomography^{†}^{†}thanks: The research of the authors was supported by AFOSR grant FA95500910091.
Abstract
In this note, we examine the problem of identifying the interaction geometry among a known number of agents, adopting a consensustype 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 steeredandobserved 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 degreebased 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 nondestructive 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 consensustype coordination algorithms. Consensustype algorithms have recently been employed for analysis and synthesis of a host of distributed protocols and control strategies in multiagent 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 informationexchange 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 inputoutput 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 degreebased graph reconstruction. The implicit contribution of our analysis is its ramifications for exact identification from boundary nodes for networks that have an embedded consensustype algorithms for their operation, including formation flying, distributed estimation, and mobile robotics.
Our notation and terminology are standard.^{1}^{1}1The 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 twoelement 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 semidefinite 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
(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 inputoutput linear timeinvariant system,
(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 inputoutput 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 graphbased coordination algorithms, namely, the feasibility of identifying the spectral and structural properties of the underlying network via the data facilitated by the inputoutput 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 inputoutput 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.
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,
(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. ∎
Iia 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,^{2}^{2}2The 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
(4) 
where , , , , , and .^{3}^{3}3The notation for a square matrix refers to its matrix exponential. In fact, the system identification process leads to a realization of the model
(5) 
where is the realization of in (4). The estimated system (5), on the other hand, is equivalent to the continuoustime system
(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).^{4}^{4}4This follows from Lemma II.1 since will appear as zero in the corresponding entries of . For example in the identification procedure called Iterative PredictionError 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 meansquare 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 inputoutput observations. The Hankel matrix, which can be constructed from the gathered inputoutput 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 continuoustime 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 steeredandobserved 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 inputoutput matrices, reduce the search space for the underlying interaction geometry? In this note, we explore this question using techniques based on integer partitioning and degreebased 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.
Iiia Graph Characterization via the Characteristic Polynomial
Recall that via a system identification method, the characteristic equation of the system (2) can be found as
(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 wellknown 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.
IiiB 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 degreebased 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
(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).
IiiC 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
(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
(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,
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.
IiiD Degree Based Graph Construction Algorithms and Complexity Analysis
Before describing the algorithm, we need to provide a few definitions.^{5}^{5}5The 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
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 nonincreasing 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;^{6}^{6}6In the case where , entries of are known. define these set of edges as being “predetermined” 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 .
Given a graphical sequence .

Define the rightmost adjacency set containing the largest index nodes different from . Let us also define as the forbidden neighbors of node . Note that originally might contain some nodes according the forbidden set . Create the set and for node connect node to (this never breaks graphicality). Set . Define the new sequence . Let .

Connect another edge of to . Run the graphicality test in Theorem III.3.

If this test fails, set ; repeat I.1.

If the test passes, keep (save) the connection, add the node to the forbidden set and update the degree sequence , set , and if has edges left, repeat from I.1.


Create the set of all adjacency sets of node that are colexicographically smaller than and preserve graphicality, i.e.,
where the operation means a colexicographic order between two sets.

For every create all graphs from the corresponding graph realization of using this algorithm, where is the sequence reduced by .
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:

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.

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

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 degreebased 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 inputoutput 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
And also
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 consensustype 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 steeredandobserved 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, “3d 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 MultiAgent Networks. Princeton University Press, Forthcoming.
 [8] M. NabiAbdolyousefi 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: PrenticeHall: 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. 13, 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, “Degreebased 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 havelhakimi type algorithm to realize graphical degree sequences of directed graphs,” arXiv:0905.4913v2, 2010.