A Passivity-Based Network Identification Algorithm with Minimal Time Complexity

A Passivity-Based Network Identification Algorithm with Minimal Time Complexity

Miel Sharf and Daniel Zelazo M. Sharf and D. Zelazo are with the Faculty of Aerospace Engineering, Israel Institute of Technology, Haifa, Israel. msharf@tx.technion.ac.il, dzelazo@technion.ac.il.
Abstract

The theory of network identification, namely identifying the (weighted) interaction topology among a known number of agents, has been widely developed for linear agents over recent years. However, the theory for nonlinear agents is far less developed, and non-applicable to large systems due to long running times. We use the notion of maximal equilibrium-independent passivity (MEIP) and network optimization theory to present a network identification method for nonlinear agents. We do so by first designing a sub-cubic time algorithm for LTI agents, and then augment it by linearization to achieve a sub-cubic time algorithm for network reconstruction for nonlinear agents and controllers. Lastly, we study the problem of network reconstruction from a complexity theory standpoint, showing that the presented algorithms are in fact optimal in terms of time complexity. We provide examples of reconstructing large-scale networks, including a network of first-order linear agents, and a non-linear neural network model.

I Introduction

Multi-agent systems have been in the pinnacle of research for the last few years, both for their variety of applications and their deep theoretical framework. They have been applied in a wide range of domains, including robotics rendezvous, formation control, flocking, distributed estimation, and social networks [1, 2, 3]. One of the most important aspects in multi-agent systems, both in theory and in practice, is the information-exchange layer, governing the interaction of agents with one another. Identifying the underlying network of a multi-agent system from measurements is of great importance in many applications. Examples include data mining and privacy in social networks [4], estimating brain connectivity using EEG in neuroscience [5], and estimating influence between currencies using a history of their exchange rates in finance [6]. Other fields with similar problems include systems biology [7, 8], communication networks [9, 10], ecology [11, 12] and physics [13, 14].

The problem of network identification has been widely studied for linear agents. Seminal works dealing with network identification include [15] in which sparse enough topologies can be identified from a small number of observations, and [16, 17], providing exact reconstruction for tree-like graphs. Other important works include [18], presenting a sieve method for solving the network identification problems for consensus-seeking networks, and [19], using a node-knockout method. More recent methods include spectral methods [20] and auto-regressive models [21]. However, a theory for network identification for interacting nonlinear agents is far less developed. We aim to provide in this work a network identification scheme for a wide range of systems, including nonlinear ones. Our approach relies on a concept widespread in multi-agent systems, namely passivity theory.

Passivity theory is a cornerstone of the theoretical framework of networks of dynamical systems [22]. Its main upshot is that it allows for the analysis of multi-agent systems to be decoupled into two separate layers, the dynamic system layer and the information exchange layer. Passivity theory was first applied to network systems by Arcak in [23]. Since then, many variations and extensions of passivity have been applied in different aspects of multi-agent systems. For instance, the related concepts of incremental passivity or relaxed co-coercivity have been used to study various synchronization problems [24, 25], and more general frameworks including Port-Hamiltonian systems on graphs [26].

One prominent variant is maximal equilibrium-independent passivity (MEIP), which was applied in [27] in order to reinterpret the analysis problem of a multi-agent system as a network optimization problem. Network optimization is a branch of optimization theory dealing with optimization of functions defined over graphs [28]. The main result of [27] showed that the asymptotic behavior of these networked systems is (inverse) optimal with respect to a family of network optimization problems. In fact, the steady-state input-output signals of both the dynamical systems and the controllers comprising the networked system can be associated to the optimization variables of either an optimal flow or an optimal potential problem; these are the two canonical dual network optimization problems described in [28]. The results of [27] were used in [29, 30] in order to solve the synthesis problem for multi-agent systems. Furthermore, the results of [27, 29] were applied in [31] to solve the network reconstruction problem for nonlinear agents. However, the presented reconstruction algorithm had super-exponential time complexity, prohibiting its application on real-life systems having more than a handful of agents.

We aim to use this network optimization framework to provide an optimal network identification scheme for multi-agent systems. We do so by injecting constant exogenous inputs, and tracking the output of the agents. By appropriately designing the exogenous inputs, we are able to reconstruct the underlying graph. The key idea in the proof is that the steady-state outputs are solutions to network optimization problems and they are one-to-one dependent on the exogenous input. This dependency can be linearized. and the associated matrix can be found by looking at a finite number of inputs and outputs. Our contributions are stated as follows:

• We present an exact sub-cubic algorithm for network reconstruction for LTI systems with MEIP agents and controllers.

• We augment the algorithm to get an approximate sub-cubic algorithm for network reconstruction for general MEIP agents and controllers.

• We explore the complexity theory behind network reconstruction algorithms, and prove that the algorithms we presented are optimal in terms of time complexity.

The rest of the paper is organized as follows. Section II surveys the relevant parts of the network optimization framework. Section III presents the problem formulation. Section IV presents the network reconstruction algorithm for LTI agents and controllers. Section V modifies the previous algorithm so it will be applicable for general MEIP agents and controllers. Section VI presents complexity bounds on general algorithms for network reconstruction, and shows that the presented algorithms are optimal. Lastly, we present a case study simulating the network reconstruction algorithms discussed in two cases - a network of LTI agents and controllers, and a nonlinear neural network.

Notations

We use basic notions from algebraic graph theory [32]. An undirected graph consists of a finite set of vertices and edges . We denote by the edge that has ends and in . For each edge , we pick an arbitrary orientation and denote . The incidence matrix of , denoted , is defined such that for edge , , , and for . Furthermore, a weighted graph is a pair where is a graph and is a collection of real numbers assigned to the edges. We say that is a positively-weighted graph if for all edges . We also use basic notations from linear algebra. We denote the standard basis vectors in by . Moreover, for a linear map between vector spaces, we denote the kernel of by , and the image of by . Furthermore, if is a subspace of an inner-product space (e.g., ), we denote the orthogonal complement of by . The notation () means the matrix is positive semi-definite (positive definite). Moreover, if , we denote the minimal eigenvalue of by . We shall also employ basic notations from complexity theory. For two functions , we say that if there exists a constant and such that for all integers . We also write if there exists a constant and such that for , i.e. . Lastly, we say that if and .

Ii Preliminaries

The role of network optimization theory in cooperative control was introduced in [27], and was used in [30] to solve the synthesis problem for multi-agent systems. In this section, we provide an overview of the main results from these works.

Ii-a The Closed-Loop and Steady-States

Consider a collection of agents interacting over a network . Assign to each node (the agents) and each edge (the controllers) the dynamical systems,

 Σi:{˙xi=fi(xi,ui)yi=hi(xi,ui),Πe:{˙ηe=ϕe(ηe,ζe)μe=ψe(ηe,ζe). (1)

We consider stacked vectors of the form and similarly for and and the operators and . The network system is diffusively coupled with the control input to each system described by , and the controller input described by . The closed-loop system is denoted by , and is is illustrated in Fig. 1.

Of interest for these systems are the steady-state solutions, if they exist, of the closed-loop. Suppose that is a steady-state of the system. Then is a steady-state input-output pair of the -th agent, and is a steady-state pair of the -th edge. This motivates the following definition, originally introduced in [27].

Definition 1.

The steady-state input-output relation of a dynamical system is the collection of all steady-state input-output pairs of the system. Given a steady-state input and a steady-state , we define

 k(u)={y: (u,y)∈k} and k−1(y)={u: (u,y)∈k}.

Let be the steady-state input-output relation for the -th agent, be the steady-state input-output relation for the -th controller, and be their stacked versions. Then, the network interconnection shown in Fig.1 imposes on the closed-loop steady-states that , , , and . As shown in [29], this is equivalent to being a steady-state for the system if and only if

 0∈k−1(y)+EGγ(ETGy).

The above expression summarizes both the algebraic and dynamic constraints that must be satisfied by the network system to achieve a steady-state solution.

Ii-B EIP Systems, MEIP Systems and Convergence of the Closed-Loop

Convergence of the system can be guaranteed under a passivity assumption on the agent and controller dynamics [27].

Definition 2 (Maximal Equilibrium Independent Passivity [27]).

Consider the dynamical system of the form

 Υ:{˙x=f(x,u)y=h(x,u), (2)

with steady-state input-output relation . The system is said to be (output-strictly) maximal equilibrium independent passive (MEIP) if the following conditions hold:

1. The system is (output-strictly) passive with respect to any steady state pair .

2. The relation is maximally monotone. That is, if then either , or , and is not contained in any larger monotone relation [33].

Such systems include simple integrators, gradient systems, Hamiltonian systems on graphs, and others (see [27, 30] for more examples). We remark that the monotonicity requirement is used to prove existence of a closed-loop steady-state, see [27] or [30] for more details.

Theorem 1 ([27, 29]).

Consider the closed-loop system . Assume that the agents are MEIP, and that the agents are output-strictly MEIP. Then the signals of the closed-loop system converge to some steady-state values satisfying .

Ii-C The Complexity of Matrix Multiplication

The time complexity of matrix multiplication is one of the most fundamental questions in complexity theory [34, 35]. The schoolbook matrix multiplication algorithm solves the problem of multiplying two matrices in time. For many years, it was believed that no faster algorithms exist.

That changed in the late 1960s when Strassen released his seminal work on matrix multiplication [36]. In this work, he exhibited a matrix multiplication algorithm that uses a divide-and-conquer method, splitting each of the two matrices to four different parts. Then, instead of multiplying these matrices blockwise, the algorithm computes seven new products, and then uses matrix addition to compute the product instead. This simple algebraic trick gives a lower time complexity, namely . When used in practice, this algorithm is a little less numerically stable then the classical algorithm, but works faster when , allowing its implementation for large matrices.

Over the next few years, algorithms with better asymptotic time complexity were found using more complex divide-and-conquer techniques. The current state-of-the-art algorithm is due to Le Gall [37], which is heavily-based on the Coppersmith-Winograd algorithm [38]. Its time complexity is about . However, the constant in front of is extremely large, namely it’s worse than even the schoolbook matrix multiplication algorithm on matrices that can be manipulated using modern-day hardware [39].

The essential time complexity of matrix multiplication is usually denoted , where poly-logarithmic terms are neglected [35]. It is widely believed that , namely that matrices can be multiplied in about time for some polynomial [35]. The current lower bound is due to Raz [40], which proved that in a specific computational model, matrix multiplication takes at least time.

It is widely known that matrix inversion and matrix multiplication have the same time complexity [34, 35]. We will take special interest in inversion of positive-definite matrices. We suspect that the time complexity of this restricted problem is the same as the time complexity of general matrix inversion, but we did not manage to find any meaningful results in the literature in this context. For the rest of this paper, we denote the time complexity of the chosen algorithm for inverting positive definite matrices by , similarly neglecting poly-logarithmic terms. This allows us to distinguish between real-world applications (in which we use the classical algorithm or Strassen’s algorithm) to theoretical complexity bounds (in which we can use the Coppersmith-Winograd algorithm). Moreover, our results still hold even if inversion of positive definite matrices turns out to be easier than inversion of general matrices. However, it should be noted that the inequality holds. Indeed, as any general matrix inversion algorithm can also be applied to positive-definite matrices. Moreover, as reading all of the input’s entries requires time.

Iii Motivation and Problem Formulation

The problem of network identification we aim to solve can be stated as follows. Given a multi-agent system , determine the underlying graph structure from the network measurements and an appropriately designed exogenous input . Many works on network identification consider networks of consensus-seeking agents [19, 18],

 ˙xi=∑{i,j}∈Eνij(xj−xi)+Biwi, (3)

where is the controlled exogenous input for the -th agent, and are the coupling coefficients. We consider a more general case of (possibly nonlinear) agents interacting over a modified protocol,

 ˙xi=fi(xi)+∑{i,j}∈Eνijgij(hj(xj)−hi(xi))+Biwi, (4)

where , and are smooth functions.111The functions are defined for all pairs, even those absent from the underlying graph. It is often assumed in multi-agent systems that each agent knows to run a given protocol (i.e., consensus). Examples of systems governed by (4), for appropriate choice of functions , include traffic control models [41], neural networks [42], and the Kuramoto model for synchronizing oscillators [43]. We let denote the stacked versions of .

In the model (3), the standard assumption is that only certain agents can be controlled using the exogenous input (i.e., is possible), and one can observe the outputs of only certain agents. To simplify the presentation, we assume that the exogenous output can be added to all agents, and that the output of all agents can be observed. In that case, we can assume without loss of generality that . We shall also denote the matrix of the coupling coefficients as .

We note that the system (4) is a special case of the closed-loop presented in Fig. 1, where the agents and the controllers are given by

 (5)

and the network is connected using the diffusive coupling and . We would like to use the mechanisms presented in Section II to establish network identification results. We make the following assumptions on the agents and controllers, allowing us to use the framework presented in Section II. With this model, we will often write the closed-loop as .

Assumption 1.

The systems , for all , are MEIP. Furthermore, the controllers , for all , are output-strictly MEIP. Moreover, The graph is positively weighted.

Assumption 2.

The inverse of the steady-state input-output relation for each agent, , is a smooth function of . Furthermore, we assume that is a twice differentiable function of , and that the derivative for all .

Assumption 2 implies that the integral function associated with [27] is smooth and . Assumption 1 implies that is strictly monotone ascending. The stronger assumption on , namely Assumption 1, made mainly to avoid heavy technical tools.

We will also consider the special case where the agents and controllers are described by linear and time-invariant (LTI) dynamics. For such systems, the input-output relation for each agent is linear and strictly monotone, and so is the function . When is an integrator, the input-output relation is given as . In these cases, is a linear function over . In particular, for some constant . We can then define the matrix such that . Similarly, we denote , where , and .

We can now formulate two fundamental problems of network detection that we will consider.

Problem 1.

Consider the network system of the form (4) satisfying Assumptions 1 and 2 with known steady-state input-output relations for the agents and controllers. Design the control inputs so that it is possible to differentiate the network system from the network system , when .

Problem 2.

Consider the network system of the form (4) satisfying Assumptions 1 and 2 with known steady-state input-output relations for the agents and controllers, but unknown network structure and coupling coefficients . Design the control inputs such that together with the output measurements of the network, it is possible to reconstruct the graph and the coupling coefficients .

In [31], a solution for Problem 1 was given for general MEIP agents and controllers, assuming that the graph is unweighted (i.e. ). The same solution also works for weighted graphs, as the proof goes word-by-word for this case as well. The time complexity of the solution in both cases was . In the same paper, a solution for Problem 2 was also given. In the case of general MEIP agents and controllers, the time complexity of the algorithm was , while for LTI agents and controllers the time complexity was . Our goal is to improve the algorithm in the general case, but we end up improving both.

Note the framework developed in [27] requires constant signals for exogenous inputs. Thus, we will consider constant , and denote them as . As in [31], we can write an equation connecting the steady-state output to the constant exogenous input .

Proposition 1.

Suppose that the system is run with the constant exogenous input . If the agents are MEIP and the controllers are output-strictly MEIP, and is the steady-state input-output relation for the agents, then the output of the system converges to a steady-state satisfying the following equation,

 w=k−1(y)+EGNg(ETGy). (6)
Proof.

The result follows directly from Theorem 1 using ), and the network connection . We note that this is an equality and not inclusion due to Assumption 2. ∎

We note that the steady-state output depends not only on the steady-state input , but also on the incidence matrix and the weights , as seen by equation (6). We wish to find a way to use this connection to reconstruct and by running the system with exogenous inputs and measuring the corresponding steady-state outputs . We first consider LTI agents and controllers, for which the equation (6) takes a linear form.

Iv An Algorithm for LTI Agents and Controllers

Suppose that our agents and controllers are LTI, i.e.,

 (7)

for some and . In this case, the steady-state input-output relations are and , where and . Thus, (6) takes the following form,

 w=Ay+EGNBETGy. (8)

In practice, the value of is known, and the value of can be measured. Moreover, the matrices are known, but the matrices and are not. We denote the unknown matrix by , and the connecting matrix by . We note that , and that is a weighted Laplacian, meaning that we can reconstruct the graph and the coupling coefficient matrix by looking at non-diagonal entries of , or of , as is diagonal.

Proposition 2.

Suppose the matrix is known. Then the graph and the coupling coefficients can be exactly reconstructed. The graph consists of all edges such that , and the coupling coefficients are .

Proof.

Directly follows from , the fact that is a weighted Laplacian with graph and weights , and the fact that , where is a diagonal matrix. ∎

Our goal now is to be able to reconstruct the matrix only from measurements. The equation connecting the steady-state output and the constant exogenous input is , where both the vectors and are known. Thus, if we have a series of measurements and , with , such that for all , and are linearly independent, we can reconstruct the matrix . Namely, we have where is the matrix whose columns are , and is the matrix whose columns are . Thus, we can solve the reconstruction problem by running the system times with different exogenous inputs and measuring the achieved steady-states, as long as we can assure that the measured steady-state outputs will be linearly independent. We can easily enforce this by considering the properties of as a linear operator.

Proposition 3.

Let be any constant, and consider the following collection of vectors :

• If , let be any collection of linearly independent vectors. Let be the solution of (8) with .

• If , let be any set of linearly independent vectors in the space orthogonal to . Let be the solution of (8) with , where for . Define with .

Then the set consists of linearly independent vectors.

Proof.

Suppose first that . In that case, it is known that is invertible [44]. Thus, sends linearly independent sets of vectors to linearly independent sets of vectors. Now suppose that . The matrix is a weighted Laplacian, meaning that the solution to the equation is not unique in , but only in , the space orthogonal to . Moreover, the matrix preserves the space , and the restricted operator is invertible. Thus, by the same reasoning as above, the vectors are linearly independent. Moreover, they are all orthogonal to . implying that the vectors are linearly independent and completing the proof. ∎

Proposition 3 suggests an algorithm for reconstruction. In the case , we choose vectors as in the theorem, run the system with them, measure the corresponding steady-state outputs, and then do a small computation involving matrix inversion, as described in the discussion prior to the proposition. In the case , we note that , independent of the values of the weighted graph . This urges us to choose and and repeat the same procedure, this time running the system only times.

However, consider Proposition 3 in case . There, the vectors are linearly independent, but the vectors are not - the last vector is equal to zero. For reasons explained later (see Remark 6), we’ll want the vectors to be linearly independent even in the case . To remedy the problem, we instead consider the matrix in the case . Moreover, we observe that and note that if is in , then . Therefore, we choose vectors as in Proposition 3 for , and also . Defining as above, we note that , and . We will implement this scheme, in which the added term of the form will be denoted by

Remark 1.

As we said, Proposition 3 gives a reconstruction scheme - choose any (or ) linearly independent vectors in the proper space, run the system (or ) times using them as inputs, measure the steady-state outputs, and then use the discussion preceding Proposition 3 to compute the graph and the weights . Instead of doing (or ) separate experiments in which we run the system with one of the -s, we can use the global asymptotic convergence of the system to use a switching signal. We use an exogenous input whose value is changed every time the system reaches its steady-state, or -close to it. See also Remark 2 about declaring when a steady-state is reached.

We conclude this chapter with an algorithm for network reconstruction, namely Algorithm 1.

Theorem 2.

Consider a diffusively coupled system . Suppose that the agents are MEIP, and that the static nonlinearities are output-strictly MEIP. Suppose furthermore that the agent and controllers are all LTI. Then

1. Algorithm 1 outputs the correct graph and coupling coefficients.

2. The time complexity of the Algorithm 1 is .

Proof.

By Propositions 2 and 3 and Remark 6, in order to prove the first claim, it’s enough to show that the vectors chosen are linearly independent, which is obvious.

As for complexity, we go over the different parts and estimate the time complexity required. The first part, before the for-loop, takes time (just to initialize for ). The first for-loop takes time as well, as each iteration takes time just to store in memory. The computation between the two for-loops takes time. Lastly, the last for-loop also takes time. Since , the result is obtained. ∎

Remark 2.

Algorithm 1 (as well as Algorithm 2 presented in the next section), like the algorithm presented in [44], runs the system with fixed exogenous inputs, and then measures the steady-state output of the closed-loop system. In practice, the exact steady-state output is achieved only asymptotically, which is both unfeasible to run, and forbids the switching input scheme of Remark 1. Therefore, we must stop the system and measure the output after a finite amount of time. We are therefore interested in understanding how long to run the system to assure sufficient proximity to the true steady-state output.

There are many ways to know the desired runtime of the system, or at least some approximation of it. One can use the storage function of the closed-loop system to estimate the distance from steady-state at each point in the run, terminating when the distance from steady-state is small enough. Another solution is to stop running the system when (or ) is small enough. Other ways to determine the runtime of the system include conducting computer-based simulations, or even intuition based on the physical time constants in the agents’ dynamics.

Another method one can use is equilibrium-independent Lyapunov exponents. For LTI systems, if we have a steady-state for the agents, it corresponds to a quadratic storage function, namely of the form for some positive numbers . We can consider the closed-loop system, for which we can write

 ˙S≤−n∑i=1ρi(hi(xi(t))−hi((xss)i))2, (9)

where are the output-passivity parameters of the agents (see [27, Theorem 3.4]). The right-hand side is bounded by . In other words, we can conclude from Lyapunov’s theorem that the closed-loop system always converges to its steady-state exponentially fast with a known exponent, no matter what steady-state it has. More exactly, the storage function decays exponentially fast with exponent , meaning that by sensing the value of the storage function at the beginning of the run, we can compute a bound on how long we need to wait until we are -close to the steady-state, namely . This method can be generalized for other, nonlinear systems as well. For example, if all the agents’ measurement functions are , the same argument works. Generically speaking, inequalities bounding the measurement function from below using the storage function of the -th agent can be used to achieve certain convergence rate guarantees, that in turn allow us to establish the needed runtime of the algorithm.

Remark 3.

After running the system, we end up with matrices and want to compute . There is another way to do the same computation, which changes its time complexity from to . Indeed, we consider . If is chosen smartly, then the product can be computed in time instead of time. Indeed, this happens if is diagonal, or more generally, if is the product of at most elementary matrices (see [45] for a definition). Indeed, if this is the case, then is also a product of at most elementary matrices, and taking a product with an elementary matrix only takes time. In this case, we want to invert the matrix , but it is positive-definite. Thus we can invert it using operations instead of operations. We show below that is the product of at most elementary matrices.

Proposition 4.

The matrix constructed in Algorithm 2 is the product of no more than elementary matrices.

Proof.

This is clear for , as , so we show it for . We run a Gaussian elimination procedure on the matrix defined by the algorithm. Each row operation corresponds to a multiplication by an elementary matrix, so it suffices to show that the procedure halts after steps. We’ll show it halts after steps. Indeed, we first consider the row operations of the form for , namely add row to row . These are total row operations, leaving all first rows unaltered, and changing the last row of the matrix to . We now divide the -th row by , which is another row operation, altering the last row to . Lastly, we apply the row operations for . These operations nullify the only nonzero off-diagonal element in each row, achieving an identity matrix. Thus, by applying a total of row operations, we transformed the matrix to the identity matrix. Thus is the product of row operations, completing the proof of the theorem. ∎

Remark 4.

We compare Algorithm 1 to the reconstruction algorithm for LTI agents and controllers suggested in [44]. Both algorithms run roughly with the same time complexity, but the presented algorithm has two advantages over the one in [44]. Firstly, it does not assume that the coupling coefficients are identical, or even known. Secondly, the algorithm in [44] can sometimes be difficult to implement due to the size of the numbers involved in it. In the presented algorithm, however, we trade this problem with inverting a matrix , which is easily invertible, eliminating numerical instability. It should also be noted that the proposed algorithm is deterministic, unlike the one presented in [44].

V An Algorithm for General MEIP Agents and Controllers Using Linearization

In general, our agents and controllers might not be LTI. However, we can still try and apply the same algorithm using linearization. As we’ll see later, we will need the following technical assumption.

Assumption 3.

For each , there are at most finitely many points at which is not twice differentiable on any bounded interval.

Heading toward linearization, we first run the system with some and get . We can now linearize the equation around . If we input , then we obtain,

 w−k−1(y)=EGNg(ETGy)≈ (10) EGNg(ETGy0)+EGN∇g(ETGy0)ETGδy,

where is the steady-state output of the network, and . More precisely, we have the following result.

Proposition 5.

Suppose that the function is twice differentiable at . Then for any small enough, the equation

 EGN∇g(ETGy0)ETGδy=δw−k−1(y)+k−1(y0)+O(∥δy∥2) (11)

holds, where .

Proof.

Immediately follows from subtracting from and using Taylor expansion up to first order, where we note that the twice differentiability assumption implies that the error of the first order approximation is . ∎

As before, injecting different signals into the system and measuring the output gives vectors, . We can use (11) to estimate the value of the matrix applied on each of . As before, we can replace one of these vectors with , as we know that it lies in the kernel of the matrix. If we knew these vectors are linearly independent, we could use the same reconstruction method as in the linear case. Thus we strive to find a method in which are linearly independent.

Theorem 3.

Let be any absolutely continuous probability measure on , and suppose that we sample according to , and let be a solution to the equation . We define vectors in the following way:

• If is differentiable at and , choose for . Define as the solution to the equation and , for . Also, set .

• Otherwise, choose for . Define as the solution to the equation and .

Suppose Assumptions 1,2, and 3 hold. If is small enough, then the set is a basis for .

Before proving Theorem 3, we state and prove a lemma.

Lemma 1.

Suppose that the same assumptions as in Theorem 3 hold. Then for any and any number , the set of all such that the solution to satisfies has measure zero.

Proof.

We consider the map defined by . The relevant set is the image of under . The assumption on implies that it is continuous and piecewise smooth, hence locally Lipschitz. Thus, is absolutely continuous, sending zero-measure sets to zero-measure sets. As has measure zero, we conclude that also has measure zero, which concludes the proof. ∎

Corollary 1.

Under the same assumptions as in Theorem 3, the function is twice differentiable at with probability 1.

We can now prove Theorem 3.

Proof.

The idea of the proof is to reduce the theorem to Proposition 3 using linearization. By Corollary 1, we know that the -probability that is not twice differentiable at is zero. Thus, we can assume this scenario does not happen.

Under this assumption, we can write the following equation connecting and ,

 δwi=∇k−1(yi)δyi+EGN∇g(ETGy0)ETGδyi+O(∥δyi∥2). (12)

It follows immediately in the case , and uses in the case . Because is small, and and are twice differentiable at and , we can conclude that . Thus, we can rewrite (12) differently:

 (13)

Let us first focus on the case . The matrix is invertible [44], meaning that is linearly independent if and only if the vectors on the left-hand side of (13) are linearly independent. However, these vectors are equal to , for some vectors satisfying , making them linearly independent for small enough, meaning that is a basis (with probability 1).

As for the case in which , we note that preserves the space orthogonal to . Moreover, when restricted to that subspace, it is an invertible map. As are orthogonal to , it’s enough to show that the former are linearly independent. As the map is invertible on the space , this is the same as saying that the vectors on the left hand side of equation (13) are linearly independent. However, these vectors are of the form , which are clearly linearly independent if is small enough. Thus is a basis for (with probability 1). This concludes the proof. ∎

Remark 5.

In the proof above, we used the fact that if and are twice differentiable at and , then . In particular, the error rate in Proposition 5 is . Moreover, we note that Remark 1 still holds for nonlinear systems, so we again may use switching inputs.

We wish to conclude this section with a proper description and analysis of the algorithm. The algorithm can be read in Algorithm 2.

Note 1.

We changed the query from the original algorithm to in the augmented algorithm. This is much more robust to the inherent noise in , arising not only from numerics, but also from the negligence of the unknown quadratic term in (11).

It’s clear that this time this is an approximate algorithm, as the quadratic error term will have an effect on the coupling coefficients. However, we can still use it to reconstruct the underlying graphs well, as will be shown later by examples. We can bound the error of algorithm, and determine its time complexity.

Theorem 4.

Consider a diffusively coupled system . Suppose that the agents are MEIP, and that the static nonlinearity is output-strictly MEIP. Moreover, suppose that assumptions 1, 2 and 3 all hold. Then:

1. Let be the matrix calculated by Algorithm 2, and define

 M=∇k−1(y0)+EGN∇g(ETGy0)ETG.

Then for any , we have

 |Mij−Mij|≤O(√nκ(1+maxi,j(νijdij)λmax(G))),

with probability 1. Thus, Algorithm 2 outputs an approximation of the graph and coupling coefficients, with probability .

2. The time complexity of the Algorithm 2 is .

We first need to prove a lemma:

Lemma 2.

Let be the matrix computed by Algorithm 2. Then the operator norm of is bounded by .

Proof.

If , then the matrix is equal to , meaning that its inverse is , and the result is clear. If , however, then is equal to , where columns are given by for and . Thus, the inverse of is equal to , so it’s enough to show that the operator norm of is bounded by . Consider the Gaussian elimination procedure applied to , as described in the proof of Proposition 4. There, we used row operations to transform to . The matrix can be computed by applying the same row operations on , in the same order. This yields a closed form for .

Indeed, we follow the same process, this time on . First, we applied the row operations for . This leaves all rows but the last unaltered, and the last row becomes . Then, we apply the map , dividing the last row by . Lastly, we applied the row operations for , adding to each of the rows but the last. Thus, the matrix is the sum of two matrices. The first is , where