A PassivityBased Network Identification Algorithm with Minimal Time Complexity
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 nonapplicable to large systems due to long running times. We use the notion of maximal equilibriumindependent passivity (MEIP) and network optimization theory to present a network identification method for nonlinear agents. We do so by first designing a subcubic time algorithm for LTI agents, and then augment it by linearization to achieve a subcubic 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 largescale networks, including a network of firstorder linear agents, and a nonlinear neural network model.
I Introduction
Multiagent 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 multiagent systems, both in theory and in practice, is the informationexchange layer, governing the interaction of agents with one another. Identifying the underlying network of a multiagent 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 treelike graphs. Other important works include [18], presenting a sieve method for solving the network identification problems for consensusseeking networks, and [19], using a nodeknockout method. More recent methods include spectral methods [20] and autoregressive 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 multiagent 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 multiagent 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 multiagent systems. For instance, the related concepts of incremental passivity or relaxed cocoercivity have been used to study various synchronization problems [24, 25], and more general frameworks including PortHamiltonian systems on graphs [26].
One prominent variant is maximal equilibriumindependent passivity (MEIP), which was applied in [27] in order to reinterpret the analysis problem of a multiagent 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 steadystate inputoutput 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 multiagent 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 superexponential time complexity, prohibiting its application on reallife systems having more than a handful of agents.
We aim to use this network optimization framework to provide an optimal network identification scheme for multiagent 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 steadystate outputs are solutions to network optimization problems and they are onetoone 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 subcubic algorithm for network reconstruction for LTI systems with MEIP agents and controllers.

We augment the algorithm to get an approximate subcubic 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 positivelyweighted 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 innerproduct space (e.g., ), we denote the orthogonal complement of by . The notation () means the matrix is positive semidefinite (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 multiagent systems. In this section, we provide an overview of the main results from these works.
Iia The ClosedLoop and SteadyStates
Consider a collection of agents interacting over a network . Assign to each node (the agents) and each edge (the controllers) the dynamical systems,
(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 closedloop system is denoted by , and is is illustrated in Fig. 1.
Of interest for these systems are the steadystate solutions, if they exist, of the closedloop. Suppose that is a steadystate of the system. Then is a steadystate inputoutput pair of the th agent, and is a steadystate pair of the th edge. This motivates the following definition, originally introduced in [27].
Definition 1.
The steadystate inputoutput relation of a dynamical system is the collection of all steadystate inputoutput pairs of the system. Given a steadystate input and a steadystate , we define
Let be the steadystate inputoutput relation for the th agent, be the steadystate inputoutput relation for the th controller, and be their stacked versions. Then, the network interconnection shown in Fig.1 imposes on the closedloop steadystates that , , , and . As shown in [29], this is equivalent to being a steadystate for the system if and only if
The above expression summarizes both the algebraic and dynamic constraints that must be satisfied by the network system to achieve a steadystate solution.
IiB EIP Systems, MEIP Systems and Convergence of the ClosedLoop
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
(2) 
with steadystate inputoutput relation . The system is said to be (outputstrictly) maximal equilibrium independent passive (MEIP) if the following conditions hold:

The system is (outputstrictly) passive with respect to any steady state pair .

The relation is maximally monotone. That is, if then either , or , and is not contained in any larger monotone relation [33].
IiC 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 divideandconquer 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 divideandconquer techniques. The current stateoftheart algorithm is due to Le Gall [37], which is heavilybased on the CoppersmithWinograd 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 modernday hardware [39].
The essential time complexity of matrix multiplication is usually denoted , where polylogarithmic 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 positivedefinite 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 polylogarithmic terms. This allows us to distinguish between realworld applications (in which we use the classical algorithm or Strassen’s algorithm) to theoretical complexity bounds (in which we can use the CoppersmithWinograd 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 positivedefinite 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 multiagent system , determine the underlying graph structure from the network measurements and an appropriately designed exogenous input . Many works on network identification consider networks of consensusseeking agents [19, 18],
(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,
(4) 
where , and are smooth functions.^{1}^{1}1The functions are defined for all pairs, even those absent from the underlying graph. It is often assumed in multiagent 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 closedloop 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 closedloop as .
Assumption 1.
The systems , for all , are MEIP. Furthermore, the controllers , for all , are outputstrictly MEIP. Moreover, The graph is positively weighted.
Assumption 2.
The inverse of the steadystate inputoutput 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 timeinvariant (LTI) dynamics. For such systems, the inputoutput relation for each agent is linear and strictly monotone, and so is the function . When is an integrator, the inputoutput 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.
Problem 2.
Consider the network system of the form (4) satisfying Assumptions 1 and 2 with known steadystate inputoutput 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 wordbyword 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 steadystate 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 outputstrictly MEIP, and is the steadystate inputoutput relation for the agents, then the output of the system converges to a steadystate satisfying the following equation,
(6) 
Proof.
We note that the steadystate output depends not only on the steadystate 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 steadystate 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 steadystate inputoutput relations are and , where and . Thus, (6) takes the following form,
(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 nondiagonal 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 steadystate 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 steadystates, as long as we can assure that the measured steadystate 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 steadystate 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 steadystate 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 steadystate, or close to it. See also Remark 2 about declaring when a steadystate 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 outputstrictly MEIP. Suppose furthermore that the agent and controllers are all LTI. Then
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 forloop, takes time (just to initialize for ). The first forloop takes time as well, as each iteration takes time just to store in memory. The computation between the two forloops takes time. Lastly, the last forloop 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 steadystate output of the closedloop system. In practice, the exact steadystate 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 steadystate 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 closedloop system to estimate the distance from steadystate at each point in the run, terminating when the distance from steadystate 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 computerbased simulations, or even intuition based on the physical time constants in the agents’ dynamics.
Another method one can use is equilibriumindependent Lyapunov exponents. For LTI systems, if we have a steadystate for the agents, it corresponds to a quadratic storage function, namely of the form for some positive numbers . We can consider the closedloop system, for which we can write
(9) 
where are the outputpassivity parameters of the agents (see [27, Theorem 3.4]). The righthand side is bounded by . In other words, we can conclude from Lyapunov’s theorem that the closedloop system always converges to its steadystate exponentially fast with a known exponent, no matter what steadystate 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 steadystate, 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 positivedefinite. 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 offdiagonal 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,
(10)  
where is the steadystate 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
(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 zeromeasure sets to zeromeasure 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 ,
(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 lefthand 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.
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.
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