Intrinsic dynamics induce global symmetry in network controllability

Controlling complex networked systems to desired states is a key research goal in contemporary science. Despite recent advances in studying the impact of network topology on controllability, a comprehensive understanding of the synergistic effect of network topology and individual dynamics on controllability is still lacking. Here we offer a theoretical study with particular interest in the diversity of dynamic units characterized by different types of individual dynamics. Interestingly, we find a global symmetry accounting for the invariance of controllability with respect to exchanging the densities of any two different types of dynamic units, irrespective of the network topology. The highest controllability arises at the global symmetry point, at which different types of dynamic units are of the same density. The lowest controllability occurs when all self-loops are either completely absent or present with identical weights. These findings further improve our understanding of network controllability and have implications for devising the optimal control of complex networked systems in a wide range of fields.

A linear time-invariant system consisting of N nodes under control can be described bẏ where x and u denote state vector and input vector, respectively, B represents control matrix, describing which nodes are controlled by u, A is the state matrix representing the interactions among nodes with a ii = 0 (i = 1, 2, · · · , N ) (no self-loops). The controllability matrix of system (S1) is defined as [S1] [B, AB, A 2 B, · · · , A N −1 B]. (S2) According to the Kalman's rank condition [S1], if and only if the controllability matrix has full rank, system (S1) is fully controllable.
If we add self-loops with identical weights w s to each node in the network A, system (S1) becomeṡ where Λ = w s I N and I N is the identity matrix of order N . In order to calculate the controllability matrix of system (S3), we first need to calculate the powers of A + w s I N , expressed as We can see that (A + w s I N ) m B is the linear combination of B, AB, A 2 B, · · · , A m B.
According to Eq. (S5) and the fact that elementary transformation keeps the rank of any matrix invariant, the controllability matrix of system (S3) can be calculated as: which holds for any w s . Therefore, system (S3) has the same controllability as that of system (S1). In other words, an arbitrary network in the absence of self-loops (presence of self-loops with zero weights) are of the same controllability as that of the network full of self-loops with identical weights.

Efficient approaches for obtaining N D based on Exact-Controllability Theory
We provide details for deriving the unified formulas for N D of both sparse and dense networks with intrinsic individual dynamics of arbitrary order. Individual dynamics are represented by dynamic units that can be integrated into the network representation, which allows us to use the exact-controllability framework to compute N D efficiently.

Dynamic units with 1st-order individual dynamics
In this case, the individual dynamics are characterized by a node with self-loop, the weight of which correspond to the single coefficient in the individual dynamics. First we consider the weight of selfloops can be either zero or w s . We derive the efficient formulas for sparse and dense networks separately.
(i) Sparse networks with random link weights. Recall that the exact-controllability framework stipulates that the minimum number of driver nodes N D is determined by the maximum geometric multiplicity where λ max is the eigenvalue corresponding to max i {µ(λ i )}. Numerical experiments suggest that due to the dominance of diagonal elements to the eigenvalues of sparse matrices with random weights, with high probability λ max will take either 0 or w s . This results in efficient computation of N D : (ii) Dense networks with identical link weight w l .
We need to first explore the eigenvalues of matrix A. If A is very dense, and all links are of identical weights w l , it is of high probability to find two of many rows are of the following form: where the two nodes corresponding to the two rows are interconnected. The corresponding rows in the (S10) We see that if λ = −w l , the two rows are identical, indicating that det(−w l I N − A) = 0 and −w l is an eigenvalue of matrix (S9).If a network is very dense, the likelihood to observe two rows with in form shown in matrix (S9) is high, numerical experiments suggest that typically λ max = −w l , resulting in an efficient computation of N D : In the presence of nonzero self-loops with weights w s , following similar analysis as shown above, the expected λ max is either −w l or w s − w l , yielding an efficient computation of N D : For more than two types of self-loops, numerical experiments suggest that λ max is the dominating elements in Λ as well. In analogy with the case with two types self-loops, we can still compute N D efficiently.
(I) For sparse networks with random weights, (II) For dense networks with identical weight w l , If one type of self-loop prevails in the network, the comparison is no longer needed and λ max is the weight w max of the self-loop, reducing the above formulas to and respectively.

The symmetry of N D
We consider the symmetry of N D for structured matrix A and in the presence of multiple types of selfloops. Without loss of generality, we denote the matrix s . According to Eq. (S13), (S18) Now we exchange two types of self-loops. There are two cases: (i) the exchange does not involve w For case (i), without loss of generality, we consider the exchange between w (2) s and w (n)

s . The matrix
Φ after the exchange becomes According to Eq. (S13), Since A is a structured matrix and both Λ − w (1) are typically so-called mixed matrices [S2].
Formally, the concept of mixed matrix is defined as follows. Let K be a subfield of a field F. For is a "structured" matrix over F, i.e., its entries are either fixed zeros or algebraically independent over We denote the row set and the column set of M as R and C, respectively, and the submatrices of Q Here, according to the rank identify of mixed matrices [S2], we have which indicates N D = N D , assuring the symmetry of N D with respect to exchanging any two types of self-loops except w (1) s .
For case (ii), assume that we exchange w (1) According to Eq. (S13), It can be simply proved that rank(Λ − w (1) . Thus according to the rank identify of mixed matrices [S2], we have The proof of case (i) and case (ii) can be immediately extended to the random distribution of different self-loops in the diagonal. Based on the rank identity of mixed matrices, we can prove the symmetry of N D as well. Taken together, N D is symmetric with respect to exchanging any two types of self-loops in systems with structured matrix A. If A is not a structured matrix, we anticipate that the symmetry of N D still holds with high probability.

Multiple types of dynamic units with arbitrary order of individual dynamics
The dynamic unit is described by a with d + 1 non-zero elements, (i.e., 1, a 0 , a 1 , · · · , a d−1 ) and in general d distinct eigenvalues (except for some pathological cases with zero measure that occur when the values of a 0 , a 1 , · · · , a d−1 satisfy certain accidental constrains). The whole network system is described by a dN × dN matrix Φ.
where λ (i) d is any one of the eigenvalues of unit type i with d dimension. Note that different eigenvalues of the same type of dynamic units are typically of the same multiplicity, so any one of them can reflect their impact on N D . The prevalence of a particular dynamic unit leads to a simplified formula where λ max d is an eigenvalue of the dynamic unit with the largest fraction. The general theory for arbitrary order of individual dynamics is degenerated to that for the special case of 1st-order individual dynamics by setting d = 1 and replacing the eigenvalues λ In analogy with the 1st-order dynamics, the symmetry of N D and the highest controllability at the global symmetry point maintain for arbitrary order of individual dynamics and network structure.

Mixture of dynamic units with different orders
We explore a more general scenario of a mixture of dynamic units of different orders. In this case, the ECT is still applicable insofar as we integrate individual dynamics and network structure into the matrix form. The efficient theory can be immediately generalized to be used for quantifying n D : where η (i) d is the weight of type i self-loop for the 1st-order individual dynamics, and is the eigenvalue of type i dynamic unit for the dth-order individual dynamics, ρ d is the fraction of dth-order individual dynamics. The prevalence of a particular dynamic unit leads to a simplified formula where η max is either the weight of self-loop or the eigenvalue of the prevailing dynamic unit.

Graphical Approach
Here, we develop a graphical approach based on maximum matching to quantify the controllability of arbitrary networks with individual dynamics of any order. The cavity method developed in statistical physics has been used to study the maximum matching problem in directed networks [S3], where the unmatched nodes are nothing but the driver nodes.

For 1st-order individual dynamics
Let's denote the number of matched nodes corresponding to a maximum matching by N m (·). For a structured matrix A, according to the structural control theory, we have On the other hand, our general efficient formulas suggest that Thus we have for structured matrix A Since N m (A) can be analytically solved by the cavity method, rank(A) can then be calculated as well.
Due to the fact that matrix rank is required in our efficient theory, Eq. (S32) connects our efficient theory and the cavity method, allowing us to calculate N D by using the tool in statistical physics.
In general, the prerequisite for Eq. (S32) is the dominance of zero in the eigenvalue spectrum of matrix A. In the presence of self-loops of identical (non-zero) weights, according the efficient formula, , where w s is the weight of the dominant self-loop. Note that the faction of self-loop with w s in matrix Φ is the same as that of zero self-loop in matrix Φ − w s I N , ensuring the prevalence of zero eigenvalue in matrix Φ − w s I N . Thus the relation between matrix rank and maximum matching still holds, yielding where N m (Φ − w s I N ) can be numerically calculated by the cavity method as well. Consequently, For multiple types of self-loops integrated with sparse networks, we have For the presence of a prevailing self-loop with weight w max s , Eq. (S34) is reduced to based on the simplified formula (S15). This analysis is also valid for dense networks with random link weights, but this scenario usually leads to the trivial result of N D = 1, regardless of the configurations of self-loops.
For the special case of dense networks with identical link weights w l , the prerequisite for employing the cavity method is violated, precluding us from using it directly. This difficulty can be overcome by considering the complement graph of matrix A + w l I N . We denote the matrix whose elements are all one by J N and thus the complement graph is given by J N − A − w l I N . According to the rank inequality, we have where the equality holds if one of the ranks in the right hand side is zero. Note that rank(J N ) = 1, which is quite close to zero as compared to N if the network size is large enough. We thus approximately have The complement graph of the original network is sparse and can be related to the maximum matching as The sparsity of complement graph enables the use of Eq. (S34) for addressing N D of dense networks with identical link weights and arbitrary types of self-loops, but the weight of self-loops is changed to w l − w s : where A is dense, and J N − A except diagonal becomes sparse, ensuring the applicability of the cavity method. In the presence of a prevailing self-loop with weight w max s , Eq. (S39) is reduced to

For high-order individual dynamics
For high-order individual dynamics, the eigenvalues of dynamic units determine N D . Due to the regular configurations of dynamic units, the relationship between the maximum matching and the rank of the dN × dN network matrix is violated, so that the cavity method for the 1st-order individual dynamics cannot be extended to the high-order cases directly. The key to implementing the cavity method lies in offering a general scheme to construct an equivalent network with 1st-order individual dynamics that shows the same N D as original system. Note that since the state matrix with high-order individual dynamics must be sparse, the role of links can be regarded as the perturbation to the eigenvalues of dynamic units. We map each dynamic unit to a single node with self-loop, the weight of which is set to be an eigenvalue of the unit's state matrix. The interactions among units in the original network can be preserved by linking the correspondent single nodes. The network construction reflects the role of dynamic units in N D and the perturbation of interactions among them in a network consisting of N nodes with self-loops, offering a general mapping of networks with high-order individual dynamics into 1st-order network systems so as to use the cavity method (S34). To be concrete, for a network system with dth-order individual dynamics, N D can be calculated via where N m (·) is the maximum matching, Φ ′ N is the reduced state matrix with N nodes and λ (i) d is the one of the eigenvalues of type-i dynamic unit's state matrix. In the presence of a prevailing dynamic unit of dth order, we denote one of the eigenvalues of the dynamic unit as λ max d , Eq. (S41) is reduced to

For a mixture of dynamic units with different orders of individual dynamics
Based on the efficient theory (S28) for a mixture of individual dynamics, we have where Φ ′ N is the reduced state matrix with N nodes and η (i) d is either the weight of self-loops for the 1st-order individual dynamics or one of the eigenvalues of dth-order individual dynamics. If there exists a dynamic unit of dth order prevails in the network, Eq. (S43) is reduced to where η max d is either the self-loop of the dynamic unit if d = 1, or one of the eigenvalues of the dynamic unit's state matrix.

Calculation of the cavity method
In Supplemental Materials of Ref. [S3], the cavity method for maximum matching of directed networks is detailed. Specifically, for a directed network with in-and out-degree distributionsP (k in ) and P (k out ), respectively, the density of driver nodes is given by where z = ⟨k⟩ is the mean degree and the generating functions G(x) and G(x) are given by The quantities w 1 , w 2 , w 1 and w 2 in Eq. (S45) can be solved by the following set of self-consistent equations: where the generating functions are defined as where the terms