Enabling Controlling Complex Networks with Local Topological Information

Complex networks characterize the nature of internal/external interactions in real-world systems including social, economic, biological, ecological, and technological networks. Two issues keep as obstacles to fulfilling control of large-scale networks: structural controllability which describes the ability to guide a dynamical system from any initial state to any desired final state in finite time, with a suitable choice of inputs; and optimal control, which is a typical control approach to minimize the cost for driving the network to a predefined state with a given number of control inputs. For large complex networks without global information of network topology, both problems remain essentially open. Here we combine graph theory and control theory for tackling the two problems in one go, using only local network topology information. For the structural controllability problem, a distributed local-game matching method is proposed, where every node plays a simple Bayesian game with local information and local interactions with adjacent nodes, ensuring a suboptimal solution at a linear complexity. Starring from any structural controllability solution, a minimizing longest control path method can efficiently reach a good solution for the optimal control in large networks. Our results provide solutions for distributed complex network control and demonstrate a way to link the structural controllability and optimal control together.


Contents
i.e., their elements are either fixed zeros or independent free parameters.
Structural Controllability: A concept introduced by Lin in 1970s [1], which refers to the ability to drive a structured LTI system to any state in its entire configuration space using only certain admissible complex system manipulations. The goal of a structural controllabiity problem is typically to allocate connections between a minimum number of external control inputs and network nodes to ensure the network controllability. A system is termed as structurally controllable if it is possible to fix the independent free parameters in A, B to certain nonzero values such that system (A, B) is controllable in the usual sense, i.e., the Kalman's Controllability Rank Condition (KCRC) is satisfied. Thus, a structurally controllable system is controllable for almost all values of free parameters except for those in some proper algebraic variety in the parameter space [2] .
Neighbors : A node x j is a neighbour of x i if there exists at least one edge (either x i → x j or x j → x i ) between them.
Local Topology Information: A node x 0 's local information is defined as the input and output degrees of all its neighbors. For the implementation of LM, we define: Unmatched/Matched node: A node that is inaccessible/accessible.
Directed control path: A matching sequence of parent-child nodes which starts at an unmatched node.
Circled control path: A matching sequence of parent-child nodes which forms into an elementary circle.
Driver node: An unmatched/inaccessible node. As will be discussed, a driver node set contains the minimal set of nodes connected to external inputs for ensuring structural controllability of the network.
Control node: A node that is directly connected to an external control input, either for controllability or minimum cost control purpose.
For the implementation of MLCP, two concepts need to be defined to avoid confusion.
Independent control node: A node that is connected to a newly added external control input.
Dependent control node: A node that is connected to an existing external control input.

Related works
As illustrated in Figure 2(a) in the main paper, the main contribution of this article is the establishment of a "link" (red line) from "structural controllability" to "optimal cost control" utilizing local-game matching (LM) algorithm and minimizing longest control path (MLCP) algorithm. We uncover that local network topology information is sufficiently rich for enabling efficient control of extremely large real-life networks.
It is also observed that, for large scale relatively dense networks, a simple scheme of randomly selecting control nodes works efficiently as a suboptimal solution. In this contribution, "local-game matching (LM)", "implicit linear quadratic regulator (ILQR)" and "minimizing longest control path (MLCP)" are three essential components establishing the red link from "structural controllability" to "optimal cost control".
The Local-game Matching is a distributed, local-information based parallel algorithm for achieving the Maximum Matching in a directed network. In the past decades, a few distributed algorithms have been developed for tackling the maximum matching problem. The headstream can be traced to Hopcroft and Karp [3], who adopted a depth-first searching method to find a maximal set of augmenting paths, based on which a maximum matching can be found. Thereafter, a few more distributed/parallel algorithms have been proposed by Israeli and Itai [4], Schieber and Moran [5], and Wu and Loui [6] et al., which however request certain global topology information such as the distance between every pair of nodes or the degree distribution of the network etc [6].
Recently, further improvements have been made and a few local information based maximum matching methods have been proposed. Specifically, Hoepman [7] developed a distributed O(n)-time algorithm to find weighted matching in general graphs; Lotker [8] proposed a more general distributed algorithm for finding the maximum matching on weighted or unweighted, static or dynamic graphs; and Mansour and Vardi [9] designed a local approximation algorithm which has a high chance to solve the maximum matching problem with a low complexity of O(log 4 (n)).
The most significant difference between LM and the existing local information based algorithms is that LM requests the least amount of local information. For example, in the state-of-art algorithms proposed in [8,9], the information obtained by every node is developed from its immediate neighborhood to its radius-r(k) neighborhood after k steps of algorithms, where r(k) is a scale function of k. The greedy algorithm proposed in [7] further allows every node to discover the topology of the whole network. The LM algorithm, on the contrary, requests only limited local information throughout the whole process, hence eliminating the need of any sophisticated multi-hop information propagations.
We now briefly review some existing results related to Implicit Linear Quadratic Regulator (ILQR).
Compared with linear quadratic regulator (LQR) [10] [11], a well-known result for optimal control, we take one step further to consider the case where the input matrix of the LTI system is selectable. On the other hand, ILQR can be viewed as related to the Linear Quadratic Gaussian (LQG) problem [12,13], one of the fundamental optimal control problems. The similarity lies in that both ILQR and LQG consider the control of uncertain linear systems. In an LQG problem, the uncertainty is due to additive white Gaussian noises [14], resulting in incomplete state information. In ILQR, on the other hand, uncertainty arises since we do not know which nodes are connected to output controllers, thus having incomplete input matrix information as well as incomplete state information. It should also be pointed out that the optimization problem of ILQR formulated in this article is non-convex, which requests careful heuristic algorithm design.
There are some other related works, including those on "pinning controllability" [15], "network synchronizability" [16] [17], and consensus or agreement in multi-agent systems [18], etc. Such studies have their main focuses on whether the system can exhibit specific spatio-temporal symmetries (for instance, whether all network nodes can follow a common time-varying trajectory or converge to a common value or achieve a common goal) and the robustness of such symmetries. For example, there are algorithms for selecting "leader agents" in stochastically forced consensus networks to minimize the mean-square deviation from the consensus. The problem we address in this paper is fundamentally different from those in existing works. We study on whether the system can fully explore its state space for system control and the cost needed when applicable.
In conclusion, the fundamental difference between existing works and our study is that we exploit the networks' local topology information to achieve the optimal control objective. Essential attributes of our proposed algorithms allow the proposed methods to be truly distributed and local-information based, and hence suitable for optimal control of large scale complex networks.

Organization of this Supplementary Information
This Supplementary  • In Section 3, the formulation of the optimization problem of ILQR is discussed in Section 3.1. The convergence of "Orthonormal-constraint-based Projected Gradient Method" (OPGM) on Stiefel manifolds is proven in Section 3.2, followed by some discussions on the control node selection in Section 3.3. Finally, comparisons between PGM and OPGM are discussed in Section 3.4.
• In Section 4, a formal description of MLCP is presented in Section 4.1, followed by a simple example of algorithm implementation illustrated in Section 4.2. The performance of MLCP on synthetic and real-life networks are summarized and discussed in Section 4.3. Finally, further discussions are provided in Section 4.4 on some details of the simulations generating the results presented in Figure   4(c) in the main paper, as well as the insights we could achieve in this figure.

The List of Notations in this Supplementary Information
The out-degree of v i in the n-th iteration The in-degree of v j in the n-th iteration The out-degree of The in-degree of The mean degree T 5 The ceiling function of the average input and output degrees of these adjacent nodes for The maximum outdegree of the network in the t-th iteration T 5 Continued on next page The maximum indegree of the network in the t-th iteration The probability that node i and node j are matched in the n-th iteration T 5 The probability that a q-order augmenting The expected probability that the matchings are foremd in the t-th iteration Continued on next page Edges set of the bipartite graph, The number of iterations of the LM algorithm Neighborhood set, The set of adjacent The set of adjacent The waiting probability All The state of Node i at time t Continued on next page  The codes of the local-game matching algorithm have been already published on GitHub (available at https://github.com/PinkTwoP/Local-Game-Matching). As presented in the main paper, in the localgame matching, every node requests one of its neighbor nodes to become its parent node, and another one to become its child node, if applicable. When a node is seeking for a match, we define the number of its unmatched child (parent) nodes, i.e., the nodes that have not yet achieved a match with a parent (child), as its u-output (u-input) degree. Denote the u-input and u-output degrees of x i as N u−in (x i ) and N u−out (x i ), respectively. Also denote the minimum u-input degree of x i 's unmatched child nodes as min{N u−in (x i )}, and the minimum u-output degree of x i 's unmatched parent nodes as min{N u−out (x i )}.
The basic strategy is that, each node without a matched child (parent) node sends a child (parent) request to its neighboring child/parent node (including itself if there exists a self-loop link to the node) with the minimum u-input (u-output) degree. Figure S1 illustrates the flowchart of the LM algorithm containing the following three components: 1) Child locating: Each node without a matched child node sends a child request to its neighboring child node (including itself if there exists a self-loop link to the node) with the minimum u-input degree min{N u−in (x i )}. When there is a tie, i.e., a node has multiple unmatched child nodes with the same minimum u-input degree, the node either holds on at a probability w or randomly breaks the tie.  Figure S1: Flowchart of implementing the LM algorithm.
2) Parent locating: Each node without a matched parent node sends a parent request to the neighboring parent node (including itself if there exists a self-loop link to the node) with the minimum u-output degree min{N u−out (x i )}. When there is a tie, i.e., a node has multiple unmatched parent nodes with the same minimum u-output degree, the node either holds on at a probability w or randomly breaks the tie.
3) Matching: When there is a match of requests, i.e., two nodes receive each other's parent/child request, a parent-child match is achieved and fixed. The child node in this match removes all its links to other parent nodes, and the u-output degrees of those parent nodes are reduced by 1. The parent node in this match removes all of its links to other child nodes, and the u-input degrees of those child nodes are reduced by 1.
The iterative request-matching operations in the LM algorithm continue until no more child/parent match can be achieved.

Examples of LM in small networks
We use two simple examples to illustrate the matching process of the LM algorithm. In the first example in Figure S2, it is interesting to observe that LM manages to reach the same final answer with the minimum number of driver nodes through a few different paths. Specifically, for node x 2 , as x 5 and x 7 have the same number of input links, they have the equal chance to be selected as x 2 's child node; and for node x 6 , as x 5 and x 9 have the same number of output links, they have the equal chance to be selected as x 6 's parent node. In LM, each node has a waiting probability w to hold on when there is a tie. As some edges may be removed in the iteration process, the tie may be broken. As illustrated in Figure S2, there x Figure S2: An example of the LM method in a small network. The nodes connected by dot lines are matched nodes and x 0 is the driver node determined by LM. Figure S3: A node x i matches with another node which does not necessarily has min{N u−in (x i )} or min{N u−out (x i )} in its neighborhood in the original network. may be 4 different cases with different tie-breaking process, which however lead to the same number of driver nodes in this example. As will be discussed later in Section 2.7, different waiting processes may lead to different results, though the differences typically are not significant.
In the second example in Figure S3, we show that while a node x i always sends a child/parent request to another node with min{N u−in (x i )} or min{N u−out (x i )}, that does not necessarily mean that it will be finally matched with the node with min{N u−in (x i )} or min{N u−out (x i )} in the original network. As shown in Figure S3, each red arrow means a matched edge (i.e., an edge with its end nodes making a match) based on a bidirectional-choice rule. For example, in (a1), node 1 and node 2 form a parent-child match because node 1 requests node 2 as its child, and node 2 also requests node 1 as its parent. The whole matching process is simple and easy to understand. Figure S3(b) illustrates a different case. In this case, some nodes have self-loop links thus these nodes may send their parent/child requests to themselves.
In Local-game matching, node 3 may request node 1 as its parent with a probability of 1/3, while node 1 will request node 2 to be its child. Meanwhile nodes 2, 4 and 5 each matches with itself in the first iteration. Finally node 1 matches with node 3, a node which does not have minimum u-input degree in node 1's neighborhood in the original network (b1). Note that the number of driver nodes obtains its minimum value, i.e., N D = 1, in this example.

Theoretical Analysis of Local-game Matching
To facilitate discussions, hereafter we denote A as the adjacent matrix, V A as the vertex set, and E A as the edge set. 2) The digraph G(A, B) contains no inaccessible nodes or dilation. Remarks: Note that a node with a self-loop edge is an accessible node.  Proof. We prove by contradiction. Assume that the complex network is not structurally controlled after using the LM algorithm. Then there must exist inaccessible nodes and dilations. When there exist inaccessible dilations, from the definition of dilation, we have ∃ S ⊆ V A , |T (S)| < |S|. In other words, we have that v j lacks child and v i lacks parent. Otherwise, if v j already has a child or v i already has a parent, (v j → v i ) should have been removed according to the LM algorithm. With the link (v j → v i ) still being there, the LM algorithm can continue, as v j can send a child request to v i and v i can send a parent request to v j . Figure S4 shows a simple example. This leads to a contradiction. Hence we prove that the LM algorithm can continue untile there is no dilation left.
For inaccessible nodes, they can be controlled by assigning new external control inputs to them which makes them become driver nodes.
From above we have the conclusion that there is no inaccessible node or dilation after implementing the LM algorithm. From Lemma 1, we have that the network is structurally controlled. There are two types of control paths in a complex network after implementing the LM algorithm: directed control paths (DCPs), the origin of each of which is an unmatched node which is taken as a driver node; and circled control paths (CCPs), in which every node has a child and a parent. Figure S5 illustrates an example, from which we can see that there are only DCPs and CCPs left in the network after implementing the LM algorithm. By connecting external control inputs to the origins of DCPs, the network can be structurally controlled. Note that for the CCP, a randomly selected node needs to be connected to an existing external controller.
Definition [19]. For any directed complex network, its digraph can be represented in a corresponding bipartite graph denoted as H(A). The bipartite graph is defined as are the set of vertexes corresponding to the N columns and rows of A respectively, and edges set Figure S6 shows an example of the above definition.
We have that an augmenting path in H(A) has the properties as follows [6]: After implementing the LM algorithm, Node 1, Node 3, Node 5 and Node 7 are the origins of four directed control paths respectively. By connecting external control inputs to these four nodes (which are colored in orange), the network can be structurally controlled. For the elementary circle, a randomly selected node needs to be connected to an existing external control input.
(1) It has an odd number of edges; (2) Its origin and end are not matched; (3) Its origin and end are on the different sides of H(A); (4) Its edges are alternatively unmatched (at least one endpoint of the edge is not matched) and matched (both two endpoints of the edge are matched); (5) By reversing the matching of the edges in the augmenting path (set the unmatched edges to be matched and set the matched edges to be unmatched), the number of matchings will be increased by one.
In Figure S6, Theorem 2. Implementing the LM algorithm on digraph of a complex network is equivalent to applying it on the corresponding bipartite graph. Proof. In digraph, the LM algorithm lets each node send a parent (child) request to its neighbor node with the minimum out-degree (in-degree); such remains unchanged in the corresponding bipartite graph.
Furthermore, in digraph where a match is identified between v i and v j , the LM algorithm removes other edges of v i and v j . Such would ensure that, in the corresponding bipartite graph, a matched node will not be matched with any other nodes, which makes the matching in the LM algorithm to be equivalent to the matching in the corresponding bipartite graph.
To make it easier to understand, we describe the LM algorithm in the bipartite graph as follows: Step (1) Each node in V + A sends a child request to its adjacent node in V − A with the minimum in-degree, while each node in V − A sends a parent request to its adjacent node in V + A with the minimum out-degree; Step Step (3) If v + i and v − j form up a match in step (2), remove all the other edges connected to v + i and v − j ; this ensures that these two nodes will not send requests to any other nodes; Step (4) Repeat step (1) to step (3), until no further matches can be formed.
We define an augmenting path with q matched edges, (q+1) unmatched edges and v + i being the second node along the path as a q-order augmenting path about v + i . We also define as the probability that a q-order augmenting path about v + i is formed when node v + i sends a child request to node v i,− We have the following Theorem.
Theorem 3. Assume that 0 ≤ p un ≤ 1 (p ma = 1 − p un ) is the expected probability that a node is unmatched till the end of the LM algorithm. We have that where v − j is one of v + i 's child nodes with an input degree of k in j . The detailed definition of α n and β n is in the proof. By letting each node send a child/parent request to its adjacent nodes with the minimum input/output degrees, the LM algorithm , regardless of whether there exist ties (defined in Section 1.2) or not.
Proof. According to Theorem 2, hereafter we prove Theorem 3 in the corresponding bipartite graph.
In other words, Figure S7 shows the connection between v + i and v − j . Now we assume that there is at least one node in unmatched till the end, the two unmatched nodes, two unmatched edges and matched form an augmenting path. This happens at a probability of (1 − p to some nodes while at least one of these nodes has at least one unmatched adjacent node till the end, the two unmatched nodes, three unmatched edges and two matched edges also form an augmenting path. This happens at an expected probability of p which have no connections with the anterior part (1-order path) of the augmenting path (viz. node v + i ). Note that some nodes may have the same unmatched adjacent nodes. Each of the overlapped nodes should be counted once only. This is illustrated in Figure S8(b). In Figure S8 . On the left are the adjacent nodes of v + i and their respective incoming links, plotted in solid lines; and on the right are adjacent nodes of v − j and their respective outgoing links, plotted in dotted lines. In this case in the picture {v + i } are matched to some nodes while at least one of these nodes has at least one unmatched adjacent node till the end, those nodes in V + A with blue and green dashed lines should be ignored. Further moving forward, we may similarly denote set C 2 as . If all nodes in ∁ C2 {v j,+ y } are matched to some nodes of which all adjacent nodes are matched to some nodes which have at least one unmatched adjacent node till the end, the two unmatched nodes, four unmatched edges and three matched edges form an augmenting path. The expected probability that this happens which have no connections with the anterior part (from 1-order path to 2-order path) of the augmenting path (viz. node v + i and node v j,+ y ) and with the overlapped part being ignored.
To generalize, the augmenting path in Figure S8(a) is 1-order augmenting path about v + i and the augmenting path in Figure S8(b) is a 2-order augmenting path about v + i . Thus, if we consider the q-order augmenting path about v + i , the probability that this augmenting path exists is Therefore, based on the assumption that there is at least one node in Back to the assumption that there is at least one node in ∁ V v + i {v − j } not matched till the end, similarly the probability that an augmenting path about v − j exists can be calculated as 1 − p where β n has analogous meaning with α n . Thus the expected probability that an augmenting path is ( Now we can also calculate the expected probabilityp j, where γ n+1 and δ n+1 have the analogous meanings with α n+1 and βn + 1 respectively.
Denote v i,− x as the node with the minimum in-degree in V v + i , and its in-degree as k in min . As the LM algorithm makes v + i send a child request to v i,− x , we have that the expected probability that an is monotonically increasing with k in j , it reaches the minimum value when For v − j , similarly we have thatp j,ap is minimized as LM algorithm makes v − j send a parent request to its adjacent node with the minimal output degree.
Further consider the case where v + i has h adjacent nodes with the same minimal input degree k in does not change its form. Thus the strategy that v + i sends a child request to one of its adjacent nodes with the minimal input degree still minimizeŝ Based on Theorem 3, we can further prove that the LM algorithm allows the system to reach the equilibrium of a static Bayesian game. Proof. Based on Theorem 2, hereafter we prove Theorem 4 in the corresponding bipartite graph, i.e.
LM algorithm guides all the nodes to make its best choice with only local information.
Since the LM algorithm uses only the local information to allow all network nodes to send requests to their selected neighbor nodes simultaneously in every iteration, the structural control problem of complex network with only local information can be viewed as a static game with incomplete information, also known as the static Bayesian game [20]. Specifically, the structural control problem to which LM algorithm is applied can be translated into the game defined as (1) N ode is the set of nodes in the digraph consisting of N nodes; (2) Ω is the set of states of nature (Here "nature" can be interpreted as the consensus of all nodes. Further discussions will be provided later in this section.). Specifically, Ω = {ω 1 , ω 2 , ω 3 }, where ω 1 is that not to be on an augmenting path is most important, ω 2 is that to be matched is very important, and ω 3 is that every node knows that only local information is available to every node in the network; (4) S i is the set of strategies of Node i, which is to request for child or parent to nodes in the neighborhood of Node i based on the type t i of Node i. For (5) u i : Ω × S → R is the payoff function of Node i, i.e., the payoff of Node i to choose strategies from S i based on the nature of Ω. The payoff of Node i is only related to the Bayesian game result, and every node is only able to estimate the expected payoff. For i is matched and not on an augmenting path 0, v + i is not mathced or on an augmenting path .
For The strategy chosen by Node i is based on its payoff function u i , which is influenced by the strategy chosen by its neighbors. An equilibrium of the game G is defined to be a Nash equilibrium of the gamê For the structual control problem of complex network with finite nodes, the equivalent game equilibrium always exists [20].
Before solving the equilibrium problem, we provide some further explanations on the "nature" to make it easier to understand. We assume that in the game G, every node is in all the states of Ω. ω 1 means that every node is egoless to reduce the number of independent external control inputs required, as having an augmenting path means that there is one more matching available which could reduce one external control input. ω 2 means that every node is selfish and anxious to be matched as far as such is possible. ω 3 means that while making a decision, every node knows that limited local information is all what the other nodes can obtain.
To prove that the LM algorithm is the equilibrium, we need to show that in the game, considering the strategies of other nodes, every node can maximize its expected payoff by sending child (parent) request to its adjacent node with minimal input (output) degree. We firstly consider the probabilitŷ (3) and (4) to request for child (parent) in order to minimize the probability of forming up an augmenting path.
Therefore, no matter which type Node i is, nature ω 1 encourages Node i to choose its neighbor node with the minimal in-degree (out-degree) to request for child (parent). This preferred strategic decision can be inferred by v + i because of nature ω 3 : to maximize the probability that v + i is matched, v + i shall firstly estimate the probability each node in V v + x . If we take into consideration the situation that there are h adjacent nodes with the same minimal out-degree (viz. the tie exists) and v i,− x would randomly choose one from these h nodes to request for parent, the probability can be revised as , which means nodes with larger in-degrees are less likely to request for parent to v + i , v + i shall choose the adjacent node with the minimal in-degree to request for child to maximize the probability of making a match. Similarly, when there are h nodes with the same minimal in-degree, v − j can estimate the probability that v j,+ y requests for child to it as Thus, v − j shall choose the node with the minimal out-degree to request for parent in order to maximize the probability of making a match.
To summarize, by sending child (parent) request to the adjacent node with the minimal input (output) degree, every node can maximize the probability that it is matched and not on an augmenting path. In the game, every node wants to maximize its expected payoff, which mathematically equals the probability that it is matched and not on an augmenting path. In other words, when Node i selects an adjacent node to send a request, it firstly minimizes the probability that it is on an augmenting path, and then maximizes the probability that it is matched. The two goals can be achieved simultaneously when Node i sends a request to the adjacent node with the minimal degree. This is the equilibrium of the static Bayesian game.
Because game G is static where all nodes make strategic decisions simultaneously in each iteration, we have the probability distribution over Ω for Node i as and Therefore, the equilibrium of game G is that every node randomly chooses one from those nodes with the minimal degree among its immediate neighbors. LM algorithm guides each node to make its best choice with the limited local information. Proof. Consider a network with N nodes and in-degree/out-degree distributions p in and p out , respectively. To simplify the proof, we start with the case where each node has only a single adjacent node with the minimum in-degree or out-degree.
In the first iteration of the LM algorithm, as each node sends a request to its adjacent node with the minimal in-degree (or out-degree), the probability p ij that node i and node j are matched is decided by the probability that node i has the minimal out-degree among neighbors of node j and the probability that node j has the minimal in-degree among neighbors of node i. Thus we have that the matching probability equals where k out i and k in j are the out-degree of node i and in-degree of node j respectively.
Let the matched nodes be ignored while calculating degree distribution in the following iterations, or in other words, they may be regarded as getting removed. Denote p in,t and p out,t as the in-degree and out-degree distributions of the network in the t-th iteration, respectively. We have . Figure S9 illustrates an example. node v y has z y adjacent nodes with the same minimal out-degree. Node v x shall thus have a probability w to hold on, and a probability (1 − w) to randomly choose one from the z x adjacent nodes with the minimal in-degree to send a request. The probability that v x sends a child request to v y thus can be calculated as Applying the binomial theorem, we have Similarly, for v y we that have the probability that v y chooses v x to request for parent is Therefore, for pair (v x → v y ), the probability that (v x → v y ) is matched is For the t-th iteration, we have For the whole network, the anticipant probabilityp t of match making is i j δ i δ j|i · p t ij , where δ i is the proportion of nodes with out-degree k out i in the network, i.e., p out (k out i ), and δ j|i is the conditional probability that nodes with out-degree k out i have adjacent nodes with in-degree k in j . Hence we havê where k out,t max and k in,t max are the maximum out-degree and maximum in-degree of the network in the t-th iteration respectively, and k t max = max{k out,t max , k in,t max }.
As k out,t max i=0 k in,t max j=0 p in,t (k in,t j )p out,t (k out,t i ) = 1, denoting P ij as p in,t (k in,t j )p out,t (k out,t i ) and Q t as k out,t max · k in,t max , we have Qt i,j P ij = 1. From Hölder inequality, we have From Equation (12),p * ,t monotonically increases when k out,t max and k in,t max decrease. As the LM algorithm removes edges in every iteration , k out,t max and k in,t max monotonically decrease with an increasing value of t.
Thusp * ,t has its minimum value when t = 1. Denotep * ,1 asp * . We have thatp * is independent of the network size according to the theorem assumption, andp * is strictly positive.
Denote the number of iterations until the algorithm stops as t end and assume that, at the end of the (t end − 1)-th iteration, there are ξ node nodes left in the network (ξ node ∈ (0, N ]), which cannot make any further matches in the t end -th iteration. We have Asp m ≥p * for m = 1, 2, . . . , t end − 1, we have that and hence the expected number of iterations is capped at t end ≤ log (1−p * ) (ξ node /N ) + 1.
For convenience of discussion, hereafter we term this upper bound for the number of iterations as t end .
Note that for every node seeking a match in each iteration, at most two operations are needed, namely to send out a parent and a child requests respectively. In the first iteration,p 1 · N nodes get matched at a cost of 2 · N operations. In the second iteration, averagelyp 2 (1 −p 1 ) · N nodes get matched, after (1 −p 1 ) · 2 · N operations. In the t-th iteration (t > 1), there are t−1 m=1 (1 −p m ) · 2 · N operations. In the last iteration, there are ξ node nodes left after (2 · ξ node ) operations. Therefore, there will be a total of For simplicity, it can be derived that the number of operations has an upper bound as follows: Thus we have that the total number of operations is bounded by 2·N p * . Asp * is independent of the size of the network, the average time complexity of LM algorithm is O(N ).
From Theorem 5, we see that the LM algorithm has a linear average time complexity. For a given w, the expected number of iterations needed is independent of the network size, but related to network nodal degree distribution.

Definitions of the networks with other distributions
Let f (x) be the probability density function of the nodal degree distribution in the network. In this contribution we add to the popular ER and B-A random networks a few more cases in which network nodal degrees follow the Chi-squared distribution, the Weibull distribution, and the Gamma distribution, respectively. These distribution functions are widely used in inferential statistics.
The probability density function f (x) of the Chi-squared distribution is defined as [21] f (x, µ) = The probability density function f (x) of the Weibull distribution [22] is where k > 0 is the shape parameter and µ > 0 is the scale parameter of the distribution. The Weibull distribution is related to a number of other probability distributions. In particular, it interpolates between the exponential distribution (k = 1) and the Rayleigh distribution (k = 2). For convenience, we let k = 1 in this contribution.
The Gamma distribution [23] belongs to a family of two-parameter continuous probability distributions. Specifically, the probability density function f (x) using the shape-scale parametrization is given by for x > 0 and k, θ > 0.
Here k and θ are the shape and scale parameters respectively, and Γ(k) is the Gamma function evaluated at k.

Performance of Local-game Matching
We test the LM method on synthetic and real-life networks. The synthetic networks include the ER model [24], the BA network [25] [26], and networks generated using Chi-squared [21], Weibull [22] and Gamma [23] distributions, respectively. Topology information about all the real-life networks we tested are available from open sources (see the reference citations in Table S3). Note that both the LM and   solution. Also note that, utilizing only the local information, the LM method steadily achieves suboptimal solutions in terms of the number of driver nodes. Table S2 summarizes the calculation results produced by the two methods, and shows that the error of the LM method is typically ≪ 1%.
Figures S10(b) and S10(c) show that the two methods produce approximately the same number of driver nodes and that they also identify the nodes with approximately the same input and output degree distributions. Thus the two methods either find approximately the same set of nodes, or find two sets of nodes with approximately the same statistical properties. The sub-optimality of the LM method in these synthetic networks is thus verified.
We test the performance of the LM method in a number of real-life networks and the results are summarized in Table S3. Data of the topologies of these networks is freely available in the cited references.
Note that the number of driver nodes identified by the LM method are consistently to be close to or equal to the optimal solutions identified by the MM method.

Empirical estimation of the number of iteration steps
Simulation results verify that the number of iterations is independent of network size. Figure S11      in which the convergence speed of the algorithm is determined by how quickly a match of requests can be achieved. This can be affected by the number of parent/child nodes possessed by each node, and is therefore related to the average nodal degree but is essentially unaffected by network size. This conclusion holds in all the synthetic networks we have tested.
The LM method requests simple local operations (i.e., child/parent locating and matching operations) in each iteration. We find that the number of iterations needed is consistently low. In all the networks that we have tested, the number of iterations k turns out to be approximately linearly proportional to the average node degree µ. Figures S12 and Figure S13 show that this conclusion holds in synthetic networks. Specifically, setting the waiting probability w = 0, we repeat the simulation in 10 independently generated synthetic networks with N = 10000 and plot the average results. The "fitted line" in each figure is generated by mean square error curve fitting. As shown in Figure S12, in all the synthetic networks, the numbers of iterations needed approximate a simple linear function of the average nodal degreek = aµ + b, where both the slope a and the intercept b are approximately 1. When the waiting probability w = 0, e.g., when w = 2 3 , as shown in Figure S13, the linear relationship between the number of iterations needed and the mean nodal degree remains valid, although the values of a and b differ.
We also test the number of iterations k in 30 real-life networks listed in Table S3, where the mean nodal degree varies approximately from µ = 2 to µ = 30 with different values of the waiting probability (w = 0, w = 0.1, w = 0.3, w = 0.5, w = 2 3 , and w = 0.8). Figures S14 and S15 show that the linear relationship between k and µ basically remains valid in real-life networks. When w increases, the slope of the fitted line also increases. This is to be expected because a large value of w indicates a longer average waiting time before the decision to randomly break the tie is made. Note that the slope of the fitted line for synthetic networks with a particular w is steadily close to the slope for real-life networks with the same w; the main difference is that the linear relationship becomes more significant in synthetic networks. This is because synthetic networks are random networks while real-life networks have distinct local structures which may significantly affect the number of iterations.

Effects of the waiting probability w.
Introducing a non-zero waiting probability w often improves the performance of the LM algorithm.
Figures S16(a) and S16(b) show that this is the case because it decreases the probability that the better of two solutions is overlooked when a tie is randomly broken. The effects of the waiting probability w needs further study, however. Although a higher w value lowers the probability of overlooking the better of two solutions in random tie breaking, the cost is that a larger number of iterations is needed. As an example, we compare the algorithm performance under different values of E(n wait ) where, for a given value of w, E(n wait ) = w 1−w , as stated in Theorem 6, denotes the expected number of iterations a node needs to wait until making its decision.  Figures S16(c) and S16(d). When the performance of LM with w = 0 becomes less satisfactorily (e.g., in the five networks selected to be presented in the figure), the use of w = 0 improves performance. To keep a good balance between algorithm performance and the number of iterations, we set E(n wait ) in the range of [2,4]. Theorem 6. In the waiting operation, if a node has a waiting probability w, the expected waiting time steps is E(n wait ) = w 1−w .
Proof. Let i = n wait . We know that i is a random variable which takes its value from 0 to +∞.    Figure S16: Having a non-zero waiting probability w helps improve LM's performance. (a) A simple example illustrating that having a non-zero waiting probability w may help improve the matching performance. As shown in (a), in the first iteration, x 3 's child nodes (x 4 and x 5 ) have the same input degree, and x 4 's parent nodes (x 1 and x 3 ) have the same output degree. If there is no waiting, x 3 and x 5 may achieve a match. If x 3 waits for one iteration, however, since x 1 and x 2 will achieve a match in the first iteration and the link between x 1 and x 4 will be removed as shown in (b) , x 3 shall match with x 4 in the second iteration. Finally, the optimal solution x 0 → x 1 → x 2 → ... → x 6 can be obtained. (c) Waiting helps improve LM's performance in an ER network (N = 5000, µ = 7). (d) Waiting helps improve LM's performance in real-life networks. where Note that So finally we obtain that 3 Implicit Linear Quadratic Regulator

Problem formulation of ILQR
As mentioned in the main paper, ILQR can be formulated as a matrix optimization problem under orthonormal boundary condition, where the objective is to drive the states from any initial state x 0 ∈ R N ×1 to converge to the origin (x(t f ) = 0) during the time interval [0, t f ] at the minimum cost, where the cost function is defined as E{ [52]; the operator E[·] takes the expectation of the argument over all realizations of the random initial state. The problem of an ILQR is to determine an input matrix B ∈ R N ×M and also design control inputs u ∈ R M×1 such that the quadratic cost of controlling a directed networkẋ(t) = Ax(t) + Bu(t), x(0) = x 0 is minimum. Since the pairs B(t), u(t) and kB(t), 1 k u(t) represent the same inputs and outputs for a nonzero k, it is natural to fix the norm of matrix B. There are different ways of doing this and each way corresponds to a distinctive boundary condition. In this paper, we consider the orthonormal boundary condition, i.e., B is restricted to be an orthonormal matrix. In this case, the minimum cost control problem is formulated as a constrained optimization problem, where the objective is to drive the system from an initial state x 0 to the origin during the time interval [0, t f ] with the minimum cost: and (A, B) is controllablė  [53,54], u(t, B) is given by where W B is the Gramian matrix This is because, with Equation (23), we have By lettingt = t f − τ , we obtain Thus for any given B satisfying that (A, B) is controllable, we can always design u(t) to drive the system from an arbitrary initial state to the origin in a fixed time interval [0, t f ]. For a given t f , the control cost function E(.) is only dependent on the selection of B. Our objective is therefore to find an optimal input matrix B such that E(B) is minimized. Before presenting orthonormal-condition-based projected gradient method (OPGM) that searches such a B iteratively, we need to consider the derivative of E(t f , B) with respect to B, i.e., From Equations (23) and (24), we obtain where tr(.) denotes matrix trace operator.
By assuming that each element of initial state x 0 = [x 01 , . . . , x 0N ] T is an identical independent distributed (i.i.d) variable with zero mean and variance 1, we obtain a constrained non-convex matrix optimization problem with the input matrix B as its variables: . (A, B) is controllable, where E x 0 x T 0 = I N . Note that (A, B) is controllable requires that M ≥ N D .

Convergence Proof of OPGM
To prove the convergence of OPGM proposed in Equations (2) and (3) in the main paper, we introduce the following lemmas first.
Lemma 3. Given x ≥ 0, for every ε > 0, there exists η > 0 such that Proof. The result can be easily obtained by applying Taylor expression. Thus the proof is omitted. and and Proof. It is easy to obtain that F T (B k ) = F (B k ) andT T k =T k . From the definition of F (B k ) and ∆B k , we have Similar to Equation (31), it is easy to obtain Equations (29) and (30). Proof (which is given in Equations (2) Our gradient descent method is to minimize E(B) and N (B) simultaneously. The convergence is proven if we can establish and From Lemma 5 and the iterative process of OPGM, we have and For N (B k+1 ), from Equations (37) and (38), we have Using Taylor expansion to expand N (B) [55], we have From Lemma 6 and Lemma 7, Equation (41) can be written as as B k B T k is symmetric and F (B k ) is symmetric and negative definite. Thus, From Equations (42) and (44), we have Thus the first inequality in (36)  This means that for any ǫ > 0, ∃K ∈ N such that ∀k > K, |N (B k )| < ǫ, and this can be also written as For E(B), we have From Lemma 2 and Lemma 3, Equation (48) can be written as By using Equation (38) we have Then applying Taylor expansion to expand E(B) [55], we have From Equation (46), ∃K ∈ N, ∀k > K, we have B T k B k = I M + o(η). Then, it follows that Thus, Likewise, for and Thus, we have With Equations (52) and (55), Equation (50) can be expressed as ∃K ∈ N, ∀k > K, From Equations (47) and (56), the convergence of E(B k ) is shown for a sufficiently small η. By combining the results shown in Equation (46), E(B k ) converges to E(B * ) with B * being an orthonormal matrix.
Remarks: The update rule fromB k+1 to B k+1 in Equation (3) of the main paper is explained as follows. We define a norm function as Note that ∀B ∈ R N ×M , B T B = I M is equivalent to N (B) = 0. Let B k+1 = ρ k ·B k+1 . It is seen that there is no exact solution ρ k satisfying that N (B k+1 ) = 0 since it is generally impossible to drag a tensor of B on the boundary constraints by just automatic scaling. In the iteration, we minimize N (B k+1 ) by Then we can obtain ρ k =

Control node selection
Note that both B 0 and B * obtained by OPGM imply that each node can be connected to each of the M controllers, which may not yield a practical solution for large-size networks. A practical and interesting problem therefore is to find the minimum cost control solution under the constraint that only a small set of nodes can be directly connected to external controllers. We term this small node set as control node set and consider the case where each control node can be connected to only a unique controller for achieving the lowest connection complexity.
For this case, we propose a very simple approach to derive the control node set from B * .
where r i = j B * ij for i = 1, ..., N . Clearly, we have max{r} = 1. The importance index of node i evaluates the relative importance of this node in achieving the minimum cost control objective. We form the control node set as the first M nodes with the largest importance index values. The corresponding connection matrix, denoted as B * , can be easily constructed: if a node i is the k−th node of the control node set, we set B * i,k = 1; otherwise, B * i,k = 0.

Comparisons between PGM and OPGM
Note that the matrix optimization model in the main paper with orthornormal boundary constraint B T B = I M is revisited from the trace boundary constraint tr(B T B) = M in one of our most recent work in [52], where projected gradient method (PGM) is proposed to locate the control nodes and then determine input matrix B.
Our finding is that the orthornormal boundary condition slightly outperforms trace boundary constraint for the problem of optimal cost control of complex networks since the control nodes selected by OPGM is more concentrated on the nodes that divide an elementary dilation equally for achieving a lower energy cost. The experiment is done in a dilation topology with N = 6, M = 3. To ensure the controllability of the networks in Figure S17, node 1 must be selected as a control node. When node 1 is fixed as a control node, there are totally 7 types of control nodes selection method which ensure the controllability of the network. Among them, control node set {node 1, node 3, node 5} divides the dilation equally, and it consumes the lowest energy. By experimenting the projected gradient algorithms under different boundary constraints repeatedly, we could obtain the estimated probability for each control node set selection method. The experimental results are presented in Table S4.
We observe that OPGM tends to divide the dilation equally with a higher probability for consuming less energy. By repeating the simulation 10000 times, the probability for selecting {node 1, node 3, node 5} is around 44.63% for OPGM, which is higher than that of PGM.
In addition, if we orthogonalize B 0 to satisfy that B T 0 B 0 = I M , the performance of OPGM can be further improved as the selected control nodes are more concentrated on node set {node 1, node 3, node 5}.
This implies that although convergence of OPGM is guaranteed under arbitrary initial B 0 , having an orthogonalized initial B 0 helps improve the performance significantly.

Minimizing Longest Control Path
For the illustration of MLCP, first locate original N LM D directed control paths (DCPs) and N C circled control paths (CCPs) using LM algorithm. Denote the length of each DCP and CCP as L i for 1 ≤ i ≤ N LM D , and C j for 1 ≤ j ≤ N C , respectively. In implementing MLCP, it is assumed that all nodes have the same information that LM obtained. The steps of MLCP are presented as follows.

MLCP Algorithm
The detailed flow of MLCP algorithm is as follows: Step 1. Randomly select a CCP (j−th path, for example), and break it by removing one of its links; then imagine that it is connected to the end node of a randomly selected LCP (i−th path, for example) to form a new DCP. Whether a physical link exists between the end node of the broken CCP and the end of the randomly selected LCP does not matter, as we shall see more clearly later. Update the length of the path to L i + C j → L i ; repeat this process until there is no CCP left. Finally we shall still have N LM D DCPs.
Step 2. Determine the assignment of the m 0 new external control inputs to these newly formed DCPs.
Suppose that each path is assigned n i control inputs, obviously we have i n i = m 0 : 1+ni attains the maximum at the i-th DCP. then n i = n i + 1 end end Step 3. For each DCP, mark the κ · Li 1+κ + 1 th node as an independent control node (refer to the concept definition section) for the i−th path where κ = 0, 1, . . . , n i . Totally there are N LM D + m 0 marked nodes and each should be connected to an external control input.
Step 4. If there is no independent control node in one of the original CCPs, randomly select one node and mark it as a dependent control node, which is connected a nearest existing external input.
As can be seen from the above discussions, only a small set of nodes are determined as control nodes and to be connected to external inputs, and each control node can only be binary connected to one of the M external input, which yields a practical solution for large-size networks. For more detailed information, please refer to the following example.

Example illustration
A simple example is presented in Figure S18 to illustrate the implementation of the MLCP algorithm.
As shown in Figure S18 Thus the network is decomposed into four segments, which include blue and brown DCPs, and green and purple CCPs. (c) At the 1st step, green CCP is selected to be broken. Edge (7 → 9) is removed. A directed dashed arrow from node 15 to node 9 is added and C 1 becomes the tail part of L 2 . (d) The purple CCP is selected to be broken. Edge (11 → 12) is removed. Node 5 is assumed to be connected to Node 12. (e) At the beginning of the 2nd step, n 1 = 0 and n 2 = 0. (f) Because L1 1+n1 = 8 > L2 1+n2 = 7, the first extra external control input is connected to L 1 . Then n 1 = 1 and L1 1+n1 = 4. (g) Because L1 1+n1 = 4 < L2 1+n2 = 7, the second extra external control input is connected to L 2 . Then n 2 = 1 and L2 1+n2 = 3.5. (h) Node 1, node 5, node 8 and node 13 are set as independent control nodes which are colored in orange. Node 12 is set as a dependent control node colored in yellow. (i) Back to the origin topology of the network, external control inputs are denoted as carmine five-pointed stars.
DCPs and two CCPs. The length of the blue DCP is 5, and the length of the brown DCP is 3. The length of the green CCP is 4, and the length of the purple CCP is 3. Now we use the MLCP algorithm to locate the external control inputs. In the first step, we randomly select one CCP. In this example as shown in Figure S18(c), we choose the green CCP denoted as C 1 .
Break it by randomly removing one of its edges. In this example, edge (7 → 9) is removed. Then imagine it is connected to the end node of a randomly selected DCP. In this example, brown DCP is selected.
An directed edge is added from the end node of the brown path to the origin of broken circle (i.e., green path). The added edge is denoted as a directed dashed line. Now the green CCP becomes the tail part of the brown DCP, which increases L 2 by C 2 . Thus now L 2 = 4 + 3 = 7.
There is still one CCP left, i.e., the purple CCP as shown in Figure S18(d). Similarly as before, we break it by randomly removing one edge of it, say, edge (11 → 12), and then add a directed dashed from node 5 to node 12. Now L 1 = 5 + 3 = 8. As there is no CCP left, the MLCP algorithm goes to the second step.
At this step, we allocate the extra external control inputs. At the beginning, the blue DCP has n 1 = 0 extra external control inputs, and the brown DCP has n 2 = 0. As shown in Figure S18(e), the blue DCP has a larger Li 1+ni . Thus the first extra external control input is connected to the blue DCP. Afterwards, n 1 = 1 and n 2 is still zero, and the brown DCP now has a larger Li 1+ni as shown in Figure  S18(f). Therefore, the second external control input is connected to the brown DCP. Now n 1 = 1 and n 2 = 1, as shown in Figure S18(g).
The third step is to distribute the external control inputs as uniformly as possible in each DCP. For the blue DCP, it has two external control inputs. The first one has to be connected to the origin of the DCP to ensure the structural controllability. The second one is set on node 5 to divide DCP as evenly as possible. Because node 12, node 10 and node 11 which originally consist the purple CCP all lack independent control inputs directly setting on them, we randomly choose one node from the purple CCP (in this example, node 12 is selected which is colored in yellow as shown in Figure S18(h)) and connect the nearest external control input in the blue DCP, which is the control input set on node 5, to node 12.
As for the brown DCP, similarly as before we set the origin and node 8 as the independent control nodes. The result is shown in Figure S18(h). All the independent control nodes are colored in orange.
At last, we allocate the four external control inputs into the network as shown in Figure S18(i).
External control inputs are denoted as carmine five-pointed stars.

Performance of Minimizing Longest Control Path
MLCP is applied to a few synthetic networks including the ER network with Poisson nodal-degree distribution [24], the BA network with power-law nodal-degree distribution [25] [26], and a large number of real-life networks. To evaluate the performance of the MLCP method, we let it be compared with a method with random connections between controllers and network nodes. In order to generate a random solution that ensures network controllability, we first apply the LM algorithm to find one set of driver nodes, which is very close to the optimal solutions found by MM as shown in Table S2 and S3, Section 2. Then we randomly choose m 0 additional nodes in the network to construct the control node set. This simple method is termed as Random Allocation Method (RAM) in the main paper. For MLCP, similarly we use the LM algorithm to generate the driver node set and then use MLCP to allocate the m 0 additional control nodes. We calculate the control costs of network corresponding to MLCP and RAM solutions, respectively.
The results on both synthetic networks and real-life networks are summarized in Tables S5-S6. Data of the topologies of the real-life networks is freely available in the cited references.
From Table S5, we observe that MLCP performs comparably against OPGM though each control node is restricted to be connected to a single external input, and it steadily outperforms RAM by an average of about four orders of magnitude. From Table S6, however, the observations are rather different. As the listed networks are relatively dense networks, the control costs obtained by MLCP are averagely only about half less than those obtained by RAM. Overall, similar conclusions have held in all our extensive simulations: MLCP significantly outperforms RAM in low-degree networks while the two methods perform almost indistinguishably as the networks become dense.

Experiments and Analysis for Figure 4(c) in the Main Paper
To demonstrate that MLCP performs much better than RAM in sparse networks, while the performances of the two methods become less distinguishable as the mean degree of the network becomes higher, we illustrate on one case in an evolving network with an increasing mean nodal degree. The results are presented in Figure 4 For a sparse network with a low mean degree, the number of driver nodes has to be relatively high to ensure structural controllability, which leads to a relatively low control cost. As the mean degree increases, the number of driver nodes needed decreases and the average/maximum length of control paths gradually increases, which interestingly drives up the control cost as we set m 0 a constant value of 80 in this set of experiments. This explains the increasing control cost with the mean degree that can be observed in Figure 4(c). When the mean degree of the network is further increased, the number of driver nodes N LM D decreases insignificantly but adding edges makes the average/maximum length of control paths become shorter, which reduces the control cost. Overall, the conclusion is that either having a larger number of controllers or a higher network density helps shorten the average/maximum control path length, and consequently lowers network control cost. As to the comparisons between MLCP and RAM, we see that MLCP outperforms RAM by up to eight orders of magnitude in sparse network, while the two methods perform comparably when the mean in-/out-degree of the network is larger than 6.