# Pinning control of networks: dimensionality reduction through simultaneous
block-diagonalization of matrices

Shirin Panahi  University of New Mexico, Albuquerque, NM, US 80131  Matteo
Lodi  University of Genoa, Via Opera Pia 11A, 16154, Genova, Italy  Marco
Storace  University of Genoa, Via Opera Pia 11A, 16154, Genova, Italy
Francesco Sorrentino<EMAIL_ADDRESS>University of New Mexico, Albuquerque,
NM, US 80131

###### Abstract

In this paper, we study the network pinning control problem in the presence of
two different types of coupling: (i) node-to-node coupling among the network
nodes and (ii) input-to-node coupling from the source node to the ‘pinned
nodes’. Previous work has mainly focused on the case that (i) and (ii) are of
the same type. We decouple the stability analysis of the target synchronous
solution into subproblems of the lowest dimension by using the techniques of
simultaneous block diagonalization (SBD) of matrices. Interestingly, we obtain
two different types of blocks, driven and undriven. The overall dimension of
the driven blocks is equal to the dimension of an appropriately defined
controllable subspace, while all the remaining undriven blocks are scalar. Our
main result is a decomposition of the stability problem into four independent
sets of equations, which we call quotient controllable, quotient
uncontrollable, redundant controllable, and redundant uncontrollable. Our
analysis shows that the number and location of the pinned nodes affect the
number and the dimension of each set of equations. We also observe that in a
large variety of complex networks, stability of the target synchronous
solution is _de facto_ only determined by a single quotient controllable
block.

> In this paper we consider dynamical networks formed of coupled identical
> oscillators and use pinning control to synchronize all of the network nodes
> on a given target time evolution. The problem of stability of the entire
> network about the target solution is then studied by linearizing the network
> dynamics about the target solution. Different from previous work, we focus
> on the general case that the node-to-node coupling among the network nodes
> is different from the input-to-node coupling from the source node to the
> ‘pinned nodes’. The stability problem is then decoupled into the stability
> of independent subsystems of the lowest dimension. Our analysis is relevant
> to the analysis of several complex networks for which we see that stability
> of the target synchronous solution is often determined by the maximum
> Lyapunov exponent associated with only one subsystem that we call ‘quotient
> controllable’.

##  I Introduction

Gaining control over collective network dynamics is a hot topic in both graph
and control theory. In particular, pinning control is a feedback control
strategy largely used for imposing synchronization or consensus in complex
dynamical networks wu2008cluster ; liu2011cluster ; su2012decentralized ;
yu2013synchronization ; liu2014synchronization ; wang2015pinning ;
liu2015synchronization ; delellis2018partial ; belykh2008cluster ;
belykh2000hierarchy ; chen2022pinning ; tang2013 ; chen2014pinning ;
liu2015pinning ; orouskhani2016optimizing ; chen2018pinning . Specifically,
one or more virtual leaders (the so-called sources) are added to the network
and define its desired trajectory. Each source directly controls only a small
fraction of the network nodes (the pinned nodes), by exerting a control action
that is a function of the pinning error vector, whose $i$-th component is
given by the difference between the output of the considered source and the
output of the $i$-th node. Pinning control of time varying networks has been
studied in dariani2011effect ; feng2016cluster ; lin2021pinning ;
shi2021synchronization . Our ability to impose a desired synchronous solution
to a dynamical network can have important applications in different fields of
science, such as in physical an:2011 , social kan:2015 , multi-agent
trentelman:2013 , and biological networks russo:2009 and pinning control is a
widely adopted solution su:2013 .

Here we consider the problem of pinning control of undirected networks of
dynamical systems, for the case in which the node-to-node connectivity is
different from the connectivity exerted on the pinned nodes, which is relevant
to a variety of realistic scenarios and engineering applications. As an
example, imagine a biological network that one wants to synchronize on a
specific time evolution, by pinning some of the network nodes. However, the
type of forcing one may be able to exert on the network is likely going to be
different from the biological node-to-node interactions between the network
nodes. This motivates the study of a problem in which the node-to-node
coupling among the network nodes is different from the coupling exerted on the
pinned nodes. Similar versions of this problem have been previously
investigated in wu2009 ; barajas2018 ; vega2018 ; alanis2021 using a Lyapunov
function (V-stability) which provides a sufficient stability condition. In
this paper, we investigate stability using linearization, which provides both
necessary and sufficient conditions.

We study the stability of the target synchronous solution in an undirected
network of coupled dynamical systems, subject to the control action of one
source node. Our goal is to reduce (through proper transformations) the
stability analysis into simpler problems, which can be analyzed independently
of one another. This problem presents mathematical challenges that require the
introduction of a specific formalism; in particular, we use the techniques for
simultaneous block diagonalization of matrices (SBD) uhlig:1973 ;
maehara2010numerical ; murota2010numerical which allow the reduction of the
stability problem in problems of lower dimension. Previous work has
successfully applied this approach to both the stability of complete
synchronization Ir:So and cluster synchronization zhang2020 . However, no
characterization was provided about the dimensions of the independent sets of
equations in which the stability problem is reduced. The dimension of these
sets of equations (which are related to the diagonal blocks of some matrices)
is the lowest and we show that it may vary based on the structure of the
network and the choice of the pinned nodes. Once these sets/blocks are found,
a maximum Lyapunov exponent (MLE) can be associated to each of them in terms
of a Master Stability Function. The target synchronous solution is stable only
when the MLEs corresponding to all the sets/blocks are negative, which
provides a simple criterion for assessing stability.

In summary, by combining MSF and SBD we study the stability of the target
solution in the considered networks subject to pinning control; this leads to
a decomposition of the stability problem into four different types of
independent equations, which we call quotient controllable, quotient
uncontrollable, redundant controllable, and redundant uncontrollable. The
insight we provide on both the sizes and ‘roles’ of the blocks is the main
contribution of this paper. As stated before, though there is previous work on
the SBD reduction, to the best of our knowledge, no paper has explained the
‘reason’ for the blocks resulting from application of the SBD decomposition.

To carry out the stability analysis, we resort to two alternative
transformations: one ($T$) based on the SBD decomposition and one ($\hat{T}$)
based on the concepts of controllable subspace and equitable clusters. We are
the first ones to apply the SBD approach to the pinning control problem and in
so doing we establish a connection between the blocks in which the stability
problem is reduced and the particular choice of the network connectivity and
of the pinned nodes. In particular, we show that the sizes of the blocks
provided by a finest SBD are the same as for the blocks generated by the
transformation $\hat{T}$, which has a clear interpretation in terms of
controllability and quotient graphs. We prove that the transformation
$\hat{T}$ also provides a finest SBD. This is another important contribution
of this paper.

The rest of the paper is organized as follows. In Sec. II we introduce the
problem of pinning control of networks. In Sec. III a stability analysis is
derived for a network with pinned nodes. The solution to the stability problem
consists of two steps presented in Secs. IV and V. First, we use the SBD
method to reduce the dimension of the stability problem. Then, by using two
appropriately defined transformation matrices, we obtain the transformation
$\hat{T}$ which decouples the stability problem into four independent blocks.
We then compare the application of the SBD transformation ($T$) with the
transformation $\hat{T}$. The effects of the selection of different pinned
nodes on multiple driven blocks is studied in Sec. VI. Application of the
theory to larger complex networks is considered in Sec. VII. Finally, the
conclusions are drawn in Sec. VIII.

##  II Problem Definition: Pinning control of networks

We consider an undirected network of $N$ coupled dynamical systems, which
evolve in time according to the following equation:

$\dot{\boldsymbol{x}}_{i}(t)=\boldsymbol{F}(\boldsymbol{x}_{i}(t))+\sum_{j=1}^{N}A_{ij}[\boldsymbol{G}(\boldsymbol{x}_{j}(t))-\boldsymbol{G}(\boldsymbol{x}_{i}(t))]\quad
i=1,\cdots,N$ (1)

where $\boldsymbol{x}_{i}$ is the $m$-dimensional state of node $i$,
$\boldsymbol{F}:R^{m}\rightarrow R^{m}$ is the function that governs the
dynamics of each node when isolated, and $\boldsymbol{G}:R^{m}\rightarrow
R^{m}$ is the node-to-node coupling function. The network topology is
described by the adjacency matrix $A$. If there is a connection between the
nodes $i$ and $j$ then $A_{ij}=A_{ji}=1$ otherwise $A_{ij}=A_{ji}=0$. Eq. (1)
can be rewritten as:

$\dot{\boldsymbol{x}}_{i}(t)=\boldsymbol{F}(\boldsymbol{x}_{i}(t))+\sum_{j=1}^{N}L_{ij}\boldsymbol{G}(\boldsymbol{x}_{j}(t))\quad
i=1,\cdots,N,$ (2)

where the Laplacian matrix $L=A-D$ and $D$ is a diagonal matrix, such that
$D_{ii}=\sum_{j}A_{ij}$. This network allows a completely synchronized
solution
$\boldsymbol{x}_{1}(t)=\boldsymbol{x}_{2}(t)=\cdots=\boldsymbol{x}_{n}(t)=\boldsymbol{x}_{s}(t)$,
which obeys,

$\dot{\boldsymbol{x}}_{s}(t)=\boldsymbol{F}(\boldsymbol{x}_{s}(t)).$ (3)

A relevant question that links control theory to graph theory is: how can
control inputs be introduced to force the network state to the target
synchronous state?

Equation (3) allows for a number of different synchronous solutions determined
by its initial condition. The problem studied in pinning control is how
control inputs generated by source nodes can be designed to enforce stability
of a particular ‘target’ synchronous solution, $\textbf{x}_{t}(t)$ produced by
the initial condition $\boldsymbol{x}_{t}^{0}$,

$\dot{\boldsymbol{x}}_{t}(t)=\boldsymbol{F}(\boldsymbol{x}_{t}(t)),\quad\quad\boldsymbol{x}_{t}(0)=\boldsymbol{x}_{t}^{0}.$
(4)

In order to address this problem, we introduce control inputs
$\textbf{u}_{i}(t)$ as pinning control signals,

$\dot{\boldsymbol{x}}_{i}(t)=\boldsymbol{F}(\boldsymbol{x}_{i}(t))+\sum_{j=1}^{N}L_{ij}\boldsymbol{G}(\boldsymbol{x}_{j}(t))+\textbf{u}_{i}(t),$
(5a) $\textbf{u}_{i}(t)=\gamma
r_{i}[\boldsymbol{H}(\boldsymbol{x}_{t}(t))-\boldsymbol{H}(\boldsymbol{x}_{i}(t))],$
(5b) $\quad i=1,\cdots,N$

where the binary scalar $r_{i}=1$ ($r_{i}=0$) if node $i$ is pinned (not
pinned), and the scalar $\gamma>0$ measures the strength of the control
coupling. Here, $\boldsymbol{H}:R^{m}\rightarrow R^{m}$ is the source-to-
pinned-node coupling function. Note that the control action is only directly
active on the pinned nodes. For instance, consider the network shown in Fig.
1. The network consists of $N=11$ nodes coupled via blue edges. A source node
that provides the target trajectory is shown in red. The control input is
directly applied (red directed edges) from the source node to the pinned nodes
$1$, $2$, and $3$, which are shown in black.

######  Definition 1.

Network with inputs. A network with inputs is represented by a pair of
$N\times N$ matrices $L$ and $R$, where the Laplacian matrix $L$ describes the
network connectivity and the diagonal matrix $R$ is such that
$R_{ii}=r_{i},i=1,..,N$. We call $\mathcal{V}$ the set of the $N$ network
nodes, $\mathcal{V}_{P}$ the set of the $s$ pinned nodes, and
$\mathcal{V}_{NP}$ the set of $\tau=N-s$ non pinned nodes,
$\mathcal{V}_{P}\cup\mathcal{V}_{NP}=\mathcal{V}$ and
$\mathcal{V}_{P}\cap\mathcal{V}_{NP}=\emptyset$.

######  Definition 2.

Extended network. Given a network with inputs, its extended network is
represented by the $(N+1)\times(N+1)$ directed Laplacian matrix $\tilde{L}$,
sorrentino2007 defined below,

$\tilde{L}=\begin{pmatrix}(L-R)_{N\times N}&\boldsymbol{r}\\\
\mathbf{0}_{1\times N}&0\end{pmatrix}$ (6)

where the column vector $\boldsymbol{r}=[r_{1},r_{2},...,r_{N}]^{T}$ and
$\mathbf{0}_{1\times N}$ is the zero row-vector of dimension $N$.

Figure 1: A network of $N=11$ nodes subject to pinning control. The source
node is shown in red. The pinned nodes (i.e., the nodes that receive the
source signal) are shown in black. The remaining network nodes are in blue.
The node-to-node network connections are in blue. Red connections carry the
target solution $\boldsymbol{x}_{t}(t)$ (Eq. (4)) from the source node to the
pinned nodes.

Previous work (see e.g., sorrentino2007 ; sorrentino2007effects ) has analyzed
stability of the target synchronous solution by considering pinned nodes with
the same coupling function as the node-to-node coupling function, i.e., with
$\boldsymbol{H}=\boldsymbol{G}$. Here instead we consider a generalization
where the effect of the source node’s time evolution on the pinned nodes is
given by a different coupling function.

##  III Stability Analysis

When all the $\boldsymbol{x}_{i}(t)$ converge on the target solution
$\boldsymbol{x}_{t}(t)$, the control inputs $\textbf{u}_{i}(t)$ converge to
zero. To study the stability of the synchronous solution, a small perturbation
$\delta\boldsymbol{x}_{i}=(\boldsymbol{x}_{i}-\boldsymbol{x}_{t})$ is
considered Pe:Ca . Linearization of Eq. (5) about the target synchronous
solution
$\boldsymbol{x}_{1}(t)=\boldsymbol{x}_{2}(t)=\cdots=\boldsymbol{x}_{n}(t)=\boldsymbol{x}_{t}(t)$
yields, see also sorrentino2007 ,

$\begin{array}[]{clc}\delta\dot{\boldsymbol{x}}_{i}(t)=D\boldsymbol{F}(\boldsymbol{x}_{t}(t))\delta\boldsymbol{x}_{i}(t)+\sum_{j=1}^{N}L_{ij}D\boldsymbol{G}(\boldsymbol{x}_{t}(t))\delta\boldsymbol{x}_{j}(t)\\\
-\gamma
r_{i}D\boldsymbol{H}(\boldsymbol{x}_{t}(t))\delta\boldsymbol{x}_{i}(t),\end{array}$
(7)

$i=1,\cdots,N$, where $D$ is the Jacobian operator.

By stacking all perturbation vectors together in one vector
$\delta\boldsymbol{X}=[\delta\boldsymbol{x}_{1}^{T},\delta\boldsymbol{x}_{2}^{T},…,\delta\boldsymbol{x}_{N}^{T}]^{T}$,
Eq. (7) can be rewritten as

$\begin{array}[]{clc}\delta\dot{\boldsymbol{X}}(t)=[I_{N}\otimes
D\boldsymbol{F}(\boldsymbol{x}_{t}(t))+L\otimes
D\boldsymbol{G}(\boldsymbol{x}_{t}(t))\\\ -\gamma R\otimes
D\boldsymbol{H}(\boldsymbol{x}_{t}(t))]\delta\boldsymbol{X}(t),\end{array}$
(8)

where $\otimes$ is the Kronecker product or direct product. Our goal is to
break the $(N\times m)$-dimensional Eq. (8) into a set of independent lower-
dimensional equations, thus performing a dimensional reduction. An inherent
difficulty is due to the fact that, with the exception of a few specific cases
HYP , it is not possible to simultaneously diagonalize both matrices $R$ and
$L$. Instead, in what follows we seek a transformation that decouples the set
of Eqs. (8), by simultaneously block diagonalizing $R$ and $L$ Ir:So .

##  IV Dimensionality reduction of the stability problem through SBD

######  Definition 3.

Simultaneous Block Diagonalization of Matrices (SBD) uhlig:1973 ;
maehara2010numerical ; murota2010numerical . Given a set of $q$ square
$N$-dimensional matrices $M_{1},M_{2},\cdots,M_{q}$, an SBD transformation is
an orthogonal square matrix $T$ (transformation matrix) with dimension $N$
such that

$T^{-1}M_{k}T=\oplus_{j=1}^{l}B^{k}_{j},\quad k=1,2,\cdots,q,$ (9)

where the symbol $\oplus$ is the direct sum of matrices, $l$ denotes the
number of blocks, each block $B^{k}_{j}$ is a square matrix with dimension
$b_{j}$ and $\sum_{j=1}^{l}b_{j}=N$.

######  Definition 4.

Finest SBD. A finest SBD is an SBD for which the resulting blocks can not be
further refined uhlig:1973 . In particular this means that the size of the
blocks can not be made smaller and that the number of blocks, $l$, is largest
among all SBDs.

######  Remark 1.

A finest SBD is unique only in the number and sizes of the irreducible blocks
uhlig:1973 .

######  Remark 2.

The blocks resulting from a finest SBD are also matrices that belong to
irreducible matrix $*$-algebras according to the structure theorem for matrix
$*$-algebras.maehara2010numerical ; murota2010numerical

Application of the SBD technique to complete synchronization of networks was
first proposed in Ref. Ir:So and only recently applied to cluster
synchronization in Ref. zhang2020 .

Once we obtain $T=SBD(R,L)$, the matrices $R$ and $L$ are transformed as
follows,

$T^{-1}LT=L_{T}=\oplus_{j=1}^{l}\hat{L}_{j},$ (10a)
$T^{-1}RT=R_{T}=\oplus_{j=1}^{l}\hat{R}_{j},$ (10b)

where $L_{T}$ is the transformed matrix $L$, $R_{T}$ is the transformed matrix
$R$, and the irreducible block pairs $(\hat{L}_{j},\hat{R}_{j})$ have the same
dimensions $b_{j}$, $j=1,...,l$. One of the main contributions of this paper
is to investigate the dimensions of the block pairs
$(\hat{L}_{j},\hat{R}_{j})$, as a result of the original network topology
(through $L$) and the choice of the pinned nodes (through $R$).

###  IV.1 Finding the matrix P

The procedure described in Ref. maehara2011algorithm to find the finest SBD
transformation, $T$, requires two steps. First, finding a matrix $P$ which
commutes with each member of the set of matrices (here, $R$ and $L$), that is
$PR=RP$ and $PL=LP$, is selected. Second, the transformation matrix $T$ is
constructed to have columns corresponding to the eigenvectors of the commuting
matrix $P$.

Here, we describe the approach to compute the matrix $P$ that commutes with
both $R$ and $L$. Without loss of generality, we apply the same permutation to
both matrices $R$ and $L$ such that the matrix $R$ can be rewritten,

${R}=\begin{pmatrix}I_{s}&0\\\ 0&\mathbf{0}_{N-s}\end{pmatrix},$ (11)

where $I_{s}$ is the identity matrix of size $s$, $\mathbf{0}_{N-s}$ is the
square zero matrix of dimension $N-s$, and $s$ is the number of pinned nodes.

######  Lemma 1.

Any matrix $P$ that commutes with the matrix $R$ in (11) has the following
block-diagonal structure,

$P=\begin{pmatrix}P_{1}&0\\\ 0&P_{2}\end{pmatrix}$ (12)

where the block $P_{1}$ has dimension $s$ and the block $P_{2}$ has dimension
$N-s$.

###### Proof.

First, we consider the matrix $R$ of Eq. (11) and a generic matrix
${P}=\big{(}\begin{smallmatrix}P_{1}&P_{12}\\\
P_{21}&P_{2}\end{smallmatrix}\big{)}$ with sub-blocks having the same
dimensions as those of the matrix ${R}$. Then, from the commutation equation
$P{R}={R}P$ we obtain

$\begin{array}[]{lcl}\begin{pmatrix}P_{1}&P_{12}\\\
P_{21}&P_{2}\end{pmatrix}\begin{pmatrix}I_{s}&0\\\
0&0_{(N-s)}\end{pmatrix}=\begin{pmatrix}I_{s}&0\\\
0&\mathbf{0}_{(N-s)}\end{pmatrix}\begin{pmatrix}P_{1}&P_{12}\\\
P_{21}&P_{2}\end{pmatrix}\\\ \\\ \begin{pmatrix}P_{1}&0\\\
P_{21}&\mathbf{0}_{(N-s)}\end{pmatrix}=\begin{pmatrix}P_{1}&P_{12}\\\
0&\mathbf{0}_{(N-s)}\end{pmatrix}.\end{array}$ (13)

Therefore, to fulfill the commutation equation $PR=RP$, the sub-blocks
$P_{12}$ and $P_{21}$ must be zeros and the matrix $P$ is block diagonal with
diagonal-blocks $P_{1}$ and $P_{2}$ in which each block has the same dimension
as the diagonal-blocks of matrix $R$. ∎

######  Corollary 1.

The transformation matrix $T$ also has the same block-diagonal structure,

$T=\begin{pmatrix}T_{1}&0\\\ 0&T_{2}\end{pmatrix},$ (14)

where $T_{1}$ is the matrix of eigenvectors of $P_{1}$ and $T_{2}$ is the
matrix of eigenvectors of $P_{2}$.

######  Lemma 2.

Application of the block diagonal transformation $T$ to the matrix ${R}$,
transforms ${R}$ back to itself, i.e., $T^{-1}{R}T=R_{T}={R}$.

###### Proof.

By considering the matrix $R$ of Eq. (11) and the matrix $T$ of Eq. (14) and
by using the fact that a block-diagonal matrix can be inverted block by block,
we have

$\begin{array}[]{lcl}T^{-1}RT=\begin{pmatrix}T_{1}^{-1}&0\\\
0&T_{2}^{-1}\end{pmatrix}\begin{pmatrix}I_{s}&0\\\
0&\mathbf{0}_{(N-s)}\end{pmatrix}\begin{pmatrix}T_{1}&0\\\
0&T_{2}\end{pmatrix}=\\\ \begin{pmatrix}I_{s}&0\\\
0&\mathbf{0}_{(N-s)}\end{pmatrix}.\end{array}$ (15)

∎

The above lemma has important consequences. In fact, as we will see in Sec.
IV.2, this will lead to a decomposition of the stability problem into
equations that can be of either one of two types: driven or undriven.

######  Remark 3.

Lemma (2) shows that $R_{T}=R$. This is true independent of any permutation of
the rows and columns of the matrix $R$. In the examples that follow we will
sometimes show the matrix $R_{T}$ in a form that corresponds to permutations
of rows and columns of the matrix $R$ in Eq. (11).

###  IV.2 Driven and Undriven blocks

As stated before, the purpose of using the SBD transformation is to break the
stability problem of Eq. (8) into a set of independent equations of lowest
dimension.

Due to the block-diagonal structure of both matrices $L_{T}$ and $R_{T}$ and
to Eq. (11), they can be decoupled into the pairs $(L_{d},R_{d})$ and
$(L_{ud},0_{N-c})$,

$L_{T}=\begin{pmatrix}L_{d}&0\\\ 0&L_{ud}\end{pmatrix}\quad\quad
R_{T}=\begin{pmatrix}R_{d}&0\\\ 0&0_{(N-c)}\end{pmatrix},$ (16)

where the square matrices $L_{d},R_{d}$ have size $c\geq s$ and the square
matrix $L_{ud}$ has dimension $N-c$. We remark that in general the blocks
$L_{d}$ and $R_{d}$ are composed of smaller diagonal blocks. In particular, if
$c=s$, then $R_{d}=I_{s}$, otherwise $R_{d}$ is equal to a diagonal matrix
with $s$ entries on the main diagonal equal to $1$ and $c-s$ entries on the
diagonal equal to $0$. In what follows we refer to $(L_{d},R_{d})$ as the
driven pair and $(L_{ud},0_{N-c})$ as the undriven pair. With an abuse of
language, we also refer to $c$ as the dimension of the driven pair.

Equation (8) can be split into the following sets of equations, one driven and
one undriven,

$\begin{split}\dot{\boldsymbol{\eta}}_{d}(t)=[I_{c}\otimes
D\boldsymbol{F}(\boldsymbol{x}_{s}(t))+L_{d}\otimes
D\boldsymbol{G}(\boldsymbol{x}_{s}(t))\\\ -\gamma R_{d}\otimes
D\boldsymbol{H}(\boldsymbol{x}_{s}(t))]\boldsymbol{\eta}_{d}(t),\end{split}$
(17a) $\displaystyle\dot{\boldsymbol{\eta}}_{ud}(t)$
$\displaystyle=\left[I_{N-c}\otimes
D\boldsymbol{F}(\boldsymbol{x}_{s}(t))+L_{ud}\otimes
D\boldsymbol{G}(\boldsymbol{x}_{s}(t))\right]\boldsymbol{\eta}_{ud}(t).$ (17b)

We then see that the original $mN$-dimensional problem of Eq. (8) has been
reduced to two independent lower-dimensional equations: the driven Eq. (17a)
with dimension $mc$ and the undriven Eq. (17b) with dimension $m(N-c)$.

Next, we attempt to gain further insights into the physical meaning of the
quantity $c$.

######  Definition 5.

Dimension of the controllable subspace. Given an $N$-dimensional pair $(A,B)$,
the Kalman controllability matrix is
$K=[B,AB,A^{2}B,...,A^{N-1}B]$ogata2010modern . The dimension of the
controllable subspace is equal to the rank of the matrix $K$. The pair $(A,B)$
is completely controllable if its rank is equal to $N$.

######  Definition 6.

Transformation $\boldsymbol{T_{c}}$. Given a $N$-dimensional canonical linear,
time-invariant system given by the pair $(A,B)$, $A=A^{T},B=B^{T}$, with
dimension of the controllable subspace equal to $h$, the transformation
$T_{c}(A,B)$ splits the system into a controllable subsystem and an
uncontrollable subsystem kailath1980linear ,

$T_{c}^{-1}AT_{c}=\begin{pmatrix}A_{C}&0\\\ 0&A_{U}\end{pmatrix}\quad\quad
T_{c}^{-1}B=\begin{pmatrix}B_{C}&0\\\ 0&0\end{pmatrix},$ (18)

where the controllable subsystem is given by the $h$-dimensional pair
$(A_{C},B_{C})$ and the uncontrollable subsystem is given by the
$(N-h)$-dimensional pair $(A_{U},0)$.

Note that the transformation applied to the matrix $A$ is a similarity
transformation while the transformation applied to the matrix $B$ is not.

######  Theorem 1.

The $c$ dimension of the driven pair $(L_{d},R_{d})$ is equal to $h$ the
dimension of the controllable subspace of the pair $(L,R)$.

###### Proof.

We prove this by contradiction. First we assume $c<h$, i.e., the dimension $c$
of the driven pair $(L_{d},R_{d})$ is less than the dimension $h$ of the
controllable subspace of the pair $(L,R)$. Since the dimension of the
controllable subspace is invariant under any similarity transformation, then
the controllable subspace of the pair $(L_{T},R_{T})$ cannot be less than the
dimension of the controllable subspace of the pair $(L,R)$, which contradicts
the assumption. Then we assume $c>h$, which indicates that the pair
$(L_{d},R_{d})$ is not controllable, since the dimension of the controllable
subspace of the pair $(L_{d},R_{d})$ is less than $c$. Then, a transformation
$T_{c}$ could be applied to the pair $(L_{d},R_{d})$ that reduces it into a
(lower-dimensional) controllable pair and an uncontrollable pair. Since the
pair $(L_{d},R_{d})$ cannot be decomposed into smaller blocks as it is part of
a finest SBD, the original assumption is not satisfied. We thus conclude
$c=h$.

∎

Henceforth $c$ denotes the dimension of the controllability subspace.

######  Theorem 2.

The undriven system (17b) is composed of $N-c$ scalar equations.

###### Proof.

Since (i) the Laplacian matrix $L$ is symmetric $L=L^{T}$ and (ii) the
transformation matrix $T=SBD(L,R)$ is orthogonal $T^{-1}=(T)^{T}$, we have
that,

$\begin{array}[]{clc}(T^{-1}LT)^{T}=T^{T}L^{T}(T^{-1})^{T}=T^{-1}LT,\end{array}$
(19)

which shows that the matrix $L_{T}=T^{-1}LT$ is symmetric. It also follows
that the blocks $L_{d}$ and $L_{ud}$ are symmetric and diagonalizable. By
diagonalizing the matrix $L_{ud}$ (and taking into account the fact that the
matrix $R_{ud}=0$), the undriven equation (17b) can be decoupled into a number
of scalar equations. ∎

According to Theorems 1 and 2, the pinning control stability problem can be
decoupled into two sets of equations: $c$ driven equations and $(N-c)$
undriven equations. Also, $c$ is equal to the rank of the controllability
matrix of the pair $({L},{R})$.

Next we seek to find $P_{1}$ and $P_{2}$ such that the matrix $P$ commutes
with the matrix ${L}$.

######  Lemma 3.

The matrices $P_{1}$ and $P_{2}$ can be found by choosing a vector in the
nullspace of the $s^{2}+\tau^{2}$-dimensional matrix
$K_{1}^{T}K_{1}+K_{2}^{T}K_{2}$ defined below.

###### Proof.

By rewriting ${L}=\big{(}\begin{smallmatrix}L_{11}&L_{12}\\\
L_{21}&L_{22}\end{smallmatrix}\big{)}$, the commutation equation $P{L}={L}P$
becomes

$\begin{pmatrix}P_{1}&0\\\ 0&P_{2}\end{pmatrix}\begin{pmatrix}L_{11}&L_{12}\\\
L_{21}&L_{22}\end{pmatrix}=\begin{pmatrix}L_{11}&L_{12}\\\
L_{21}&L_{22}\end{pmatrix}\begin{pmatrix}P_{1}&0\\\ 0&P_{2}\end{pmatrix},$
(20a) $\begin{pmatrix}P_{1}L_{11}&P_{1}L_{12}\\\
P_{2}L_{21}&P_{2}L_{22}\end{pmatrix}=\begin{pmatrix}L_{11}P_{1}&L_{12}P_{2}\\\
L_{21}P_{1}&L_{22}P_{2}\end{pmatrix}.$ (20b)

Eq. (20b) corresponds to the following four equations,

$P_{1}L_{11}=L_{11}P_{1}$ (21a) $P_{2}L_{22}=L_{22}P_{2}$ (21b)
$P_{1}L_{12}=L_{12}P_{2}$ (21c) $P_{2}L_{21}=L_{21}P_{1}$ (21d)

Define the operator $\text{vec}:\mathbb{R}^{n\times m}\mapsto\mathbb{R}^{nm}$
which takes as input a matrix and returns a vector consisting of the columns
of the matrix stacked on top of each other. Note that $P_{1}$ has dimension
$s$ and $P_{2}$ has dimension $\tau=N-s$. The four equations in Eq. (21) can
be expressed as the following two linear systems of equations.

$\begin{array}[]{lcl}\left[\begin{array}[]{cc}I_{s}\otimes
L_{11}-L_{11}\otimes I_{s}&O_{s^{2},\tau^{2}}\\\
O_{\tau^{2},\tau^{2}}&I_{\tau}\otimes L_{22}-L_{22}\otimes
I_{\tau}\end{array}\right]\left[\begin{array}[]{c}\text{vec}(P_{1})\\\
\text{vec}(P_{2})\end{array}\right]\\\
=\left[\begin{array}[]{c}\bm{0}_{s^{2}}\\\
\bm{0}_{\tau^{2}}\end{array}\right]\end{array}$ (22a)
$\left[\begin{array}[]{cc}L_{12}^{T}\otimes I_{s}&-I_{\tau}\otimes L_{12}\\\
-I_{s}\otimes L_{21}&L_{21}^{T}\otimes
I_{\tau}\end{array}\right]\left[\begin{array}[]{c}\text{vec}(P_{1})\\\
\text{vec}(P_{2})\end{array}\right]=\left[\begin{array}[]{c}\bm{0}_{s\tau}\\\
\bm{0}_{s\tau}\end{array}\right]$ (22b)

By inspecting Eq. (22), the vector
$\bm{p}=\left[\begin{array}[]{cc}\text{vec}(P_{1})^{T}&\text{vec}(P_{2})^{T}\end{array}\right]^{T}$
must lie in the intersection of the nullspaces of the two matrices that appear
in Eq. (22). Let $K_{1}\in\mathbb{R}^{s^{2}+\tau^{2}\times s^{2}+\tau^{2}}$
and $K_{2}\in\mathbb{R}^{2s\tau\times s^{2}+\tau^{2}}$ be the matrices that
appear in the first and second lines of Eq. (22), respectively. In the
Supplementary Material Note 3, we show that $\mathcal{N}(K_{1})\cap
N(K_{2})=\mathcal{N}(K_{1}^{T}K_{1})\cap\mathcal{N}(K_{2}^{T}K_{2})=\mathcal{N}(K_{1}^{T}K_{1}+K_{2}^{T}K_{2})$,
where the symbol $\mathcal{N}(M)$ indicates the null subspace of the matrix
$M$. Then, we create the vector $\bm{p}=\sum_{i=1}^{n_{z}}\alpha_{i}v_{i}$
where $n_{z}$ is dimension of the nullspace of the matrix
$S=K_{1}^{T}K_{1}+K_{2}^{T}K_{2}$ and $v_{i}$ are a basis for the nullspace of
$S$, that is, eigenvectors corresponding to its zero eigenvalues. With
$\bm{p}$, form the two matrices $P_{1}$ and $P_{2}$ and construct the full
matrix $P$ that appears in Eq. (20a). Additional details of the derivation of
Eq. (22) can be found in Appendix IX.1.
∎

######  Remark 4.

The procedure outlined to find the block diagonal matrix $P$ is a
specialization of the method presented in SBD2 to simultaneously block
diagonalize two matrices when one of the matrices has the form of $R$ in Eq.
(11). Using the method of SBD2 directly on $L$ and $R$ requires finding a
vector in the nullspace of a square matrix of dimensions $(s+\tau)^{2}$ while
the method outlined in Lemma 3 requires finding a vector in the nullspace of a
square matrix of dimension $s^{2}+\tau^{2}$. As $s^{2}+\tau^{2}<(s+\tau)^{2}$
for $s,\tau\geq 1$, the method presented here is more efficient than the
method of SBD2 by exploiting the structure of $R$.

To illustrate the procedure described in this section, consider the network
with inputs with $N=5$ nodes shown in Fig. 2 (a).

Figure 2: (a) A $5$-dimensional network with inputs. Nodes are color-coded
according to the equitable clusters (see Definition 7) to which they belong.
The pinned node (node number $1$) in yellow is in a cluster, while all the
remaining nodes (node $2$ to $5$) in blue are in another cluster. The effect
of the pinning control action is graphically shown as a red arrow pointing at
node $1$. Matrices (b) and (c) show the pair $(L,R)$ corresponding to this
network and matrices (d) and (e) show the pair $(L_{T},R_{T})$ obtained by
application of the SBD transformation $T$. The pair $(L_{T},R_{T})$ contains
driven blocks (inside black boxes) and undriven blocks (inside blue boxes.)

To compute the transformation matrix $T$ for the pair $(L,R)$ shown in Fig. 2
(b,c), we first find the block diagonal matrix $P$ using Lemma 1 and Lemma 3
and obtain:

$P=\begin{pmatrix}0.6828&0&0&0&0\\\ 0&0.0036&0.7265&0.0801&-0.1274\\\
0&0.7265&-0.6427&0.5189&0.0801\\\ 0&0.0801&0.5189&-0.6427&0.7265\\\
0&-0.1274&0.0801&0.7265&0.0036\end{pmatrix}$ (23)

Now by constructing the transformation matrix $T$ from the eigenvectors of the
matrix $P$ and following Eq. (10) we obtain the pair of matrices
$(L_{T},R_{T})$ which are shown in Fig. 2 (d,e): The driven part (black blocks
in Fig. 2) has dimension $2$ because the controllability subspace of the pair
($L,R$) has dimension $c=2$.

(a) (b) (c) (d)

Figure 3: MLEs related to ($\hat{L}_{1},\hat{R}_{1}$) (blue line),
($\hat{L}_{2},\hat{R}_{2}$) (red), ($\hat{L}_{3},\hat{R}_{3}$) (yellow),
($\hat{L}_{4},\hat{R}_{4}$) (violet), for the example in Fig. 2 when the
target trajectory is either the stable chaotic attractor exhibited by a single
oscillator (a) or a 1-cycle UPO (b), or a 4-cycle UPO (c), or a 8-cycle UPO
(d).

As an example to demonstrate the use of the MLE to determine the stability of
the system, we make each of the nodes in the network in Fig. 2 a Rössler
oscillator so that $m=3$ where the state vector of each node is governed by
the function

$\boldsymbol{F}(\boldsymbol{x}_{i})=\begin{bmatrix}-y_{i}-z_{i}\\\
x_{i}+ay_{i}\\\ b+(x_{i}-c)z_{i}\\\ \end{bmatrix}$ (24)

where $a=0.1$, $b=0.1$ and $c=15$. With these parameters the oscillator
evolves within a stable chaotic attractor. The network node-to-node coupling
function is $\boldsymbol{G}(\boldsymbol{x}_{i})=[x_{i},0,0]^{T}$, whereas the
chosen input-to-node coupling function is
$\boldsymbol{H}(\boldsymbol{x}_{i})=[0,y_{i},0]^{T}$. Fig. 3 shows the MLEs
associated to Eqs. (17a) and (17b) plotted vs $\gamma$ for the network of Fig.
2.

Fig. 3(a) shows that the nodes of the network synchronize on the target
trajectory for $\gamma\in[0.9,2.3]$. In particular, it is apparent that the
driven equations (17a) depend on $\gamma$, whereas the undriven equations
(17b) are independent.

(a) (b)

Figure 4: Time plots for $\gamma=1.5$. (a) Evolution of the state variables
$x_{i}$ (solid lines, $i=1$ blue, $i=2$ red, $i=3$ yellow, $i=4$ purple, $i=5$
green) and of the first component of the target trajectory (dashed line). (b)
Evolution of $E_{i}(t)=||\boldsymbol{x}_{i}-\boldsymbol{x}_{t}||_{2}$; color
code as in panel (a).

As a harder benchmark, we assign an unstable periodic orbit (UPO) embedded in
the attractor as the target trajectory which can be selected by fine-tuning
the initial condition of Eq. (3).111We computed the UPO by using the
continuation tool MATcont dhooge:2003 ; dhooge2006matcont . According to the
theory hunt:1996 , we expect node trajectories to converge to a target
trajectory corresponding to a UPO only for a relatively high period: while
Figs. 3(b,c) show that it is not possible to have convergence to an unstable
1-cycle or 4-cycle due to the presence of a positive MLE for any $\gamma$
value, Fig. 3(d) shows that the nodes of the network synchronize on the target
trajectory (unstable 8-cycle) for $\gamma\in[1\quad 2.2]$. This is confirmed
by Fig. 4a, which shows the time evolution of the state variables $x_{i}$
(solid line) and of the corresponding component (dashed line) of the target
trajectory $\textbf{x}_{t}$ (unstable 8-cycle) for $\gamma=1.5$ and for random
initial conditions. Fig. 4b shows the time evolution of the Euclidean norm
$E_{i}$ of the difference $\boldsymbol{x}_{i}-\boldsymbol{x}_{t}$, from which
we see that $\boldsymbol{x}_{i}$ converges to $\boldsymbol{x}_{t}$ for large
time.

##  V Dimensionality reduction through an alternative transformation $\hat{T}$

In what follows we discuss two transformations, each one of which results in a
simultaneous block diagonalization (although not necessarily finest) of the
pair $(L,R)$: one is the transformation $T_{c}$ which decouples the problem
(8) into a controllable equation and an uncontrollable equation; the other one
is the transformation associated with the quotient network
sanchez2020exploiting ; golubitsky2012singularities ; golubitsky2003symmetry ,
which in turn decouples the problem (8) into two sets of equations, which will
be referred to as quotient and redundant. By applying both transformations to
the pair $(L,R)$, we obtain four independent equations: quotient controllable,
quotient uncontrollable, redundant controllable, redundant uncontrollable.
Interestingly, this decomposition into 4 independent sets of equations has the
same structure as the decomposition produced by a finest SBD.

###  V.1 Matrix $T_{c}$

From definition 6, it follows that there is a transformation $T_{c}(L,R)$,
which divides the system given by the pair $(L,R)$ in a controllable subsystem
and an uncontrollable subsystem, see also martini:2010 ,

$T_{c}^{-1}LT_{c}=\begin{pmatrix}L_{C}&0\\\ 0&L_{U}\end{pmatrix}\quad\quad
T_{c}^{-1}R=\begin{pmatrix}R_{C}&0\\\ 0&0\end{pmatrix},$ (25)

where the controllable subsystem is given by the $c$-dimensional pair
$(L_{C},R_{C})$ and the uncontrollable subsystem is given by the
$(N-c)$-dimensional pair $(L_{U},0)$.

######  Remark 5.

The first $c$ linearly-independent column vectors of the matrix $T_{c}$ can be
chosen as an orthonormal basis for the range of the Kalman controllability
matrix $K$. The remaining $N-c$ columns of $T_{c}$ are chosen such that
$T_{c}$ is a basis for $\mathbb{R}^{N}$.

######  Lemma 4.

The matrix $T_{c}^{-1}RT_{c}$ has the same block structure as the matrix
$T_{c}^{-1}R$ in Eq. (25).

###### Proof.

Due to the fact that matrix $L$ is symmetric, the transformation $T_{c}$ can
be constructed to be an orthogonal matrix $T_{c}^{-1}=T_{c}^{T}$ which can be
splitted into four sub-matrices as:

$T_{c}=\begin{pmatrix}T_{C1}&T_{C2}\\\ T_{C3}&T_{C4}\end{pmatrix}$ (26)

where the two square sub-matrices $T_{C1}\in\mathbb{R}^{s\times s}$ and
$T_{C4}\in\mathbb{R}^{N-s\times N-s}$ and the two rectangular sub-matrices
$T_{C2}\in\mathbb{R}^{s\times N-s}$ and $T_{C3}\in\mathbb{R}^{N-s\times s}$.
Now according to Remark 5 and by considering that (i) the matrix $R$ (11) is
symmetric and (ii) $T_{c}^{-1}R$ has the special block-diagonal form of Eq.
(25), we have

$T_{c}^{T}R=\begin{pmatrix}T_{C1}^{T}&T_{C3}^{T}\\\
T_{C2}^{T}&T_{C4}^{T}\end{pmatrix}\begin{pmatrix}R_{1}&0\\\
0&0\end{pmatrix}=\begin{pmatrix}T_{C1}^{T}R_{1}&0\\\ T_{C2}^{T}R_{1}&0\\\
\end{pmatrix}.$ (27)

We also note that based on definition 6, $T_{C2}^{T}R_{1}=0$.

From Eq. (11), we know that $R_{1}$ is equal to the identity matrix $I_{s}$,
which implies that $T_{C2}^{T}=0$. Then for the matrix $T_{c}^{T}RT_{c}$ we
have:

$T_{c}^{T}RT_{c}=\begin{pmatrix}R_{C}&0\\\
0&0\end{pmatrix}\begin{pmatrix}T_{C1}&0\\\
T_{C3}&T_{C4}\end{pmatrix}=\begin{pmatrix}R_{C}T_{C1}&0\\\ 0&0\end{pmatrix},$
(28)

which implies that the matrix $T_{c}^{T}RT_{c}$ has the same block structure
as the matrix $T_{c}^{T}R$ in Eq. (25).

∎

It should be noted that the transformation $T_{c}$ is not unique. In this
paper we build $T_{c}=[U_{c}\quad X_{N-c}]$ as follows. We define the singular
value decomposition (SVD) of $K$ as $K=U\Sigma V^{T}$. $U_{c}$ contains the
first $c$ columns of $U$ and $X_{N-c}$ contains the $N-c$ eigenvectors of $L$
that do not belong to the controllability subspace.

###  V.2 Matrix $T_{q}$

We now introduce another transformation, which we refer to as the quotient
transformation, which we show also leads to simultaneously block-diagonalizing
the pair $(L,R)$.

######  Definition 7.

Equitable cluster partition. Given the network with inputs (defined in (6))
described by the pair $({L},R)$, one can partition the set of the network
nodes $\mathcal{V}$ into subsets, i.e., equitable clusters,
$\mathcal{C}_{1},\mathcal{C}_{2},..,\mathcal{C}_{M}$,
$\cup_{i=1}^{M}\mathcal{C}_{i}=\mathcal{V}$,
$\mathcal{C}_{i}\cap\mathcal{C}_{j}=\emptyset$ for $i\neq j$, where

$\sum_{k\in\mathcal{C}_{\ell}}\tilde{L}_{ik}=\sum_{k\in\mathcal{C}_{\ell}}\tilde{L}_{jk},\quad\begin{aligned}
\forall i,j&\in\mathcal{C}_{k}\\\
\forall\mathcal{C}_{k},\mathcal{C}_{\ell}&\subset\mathcal{V}.\end{aligned}$
(29)

######  Remark 6.

Each of the equitable clusters can either contain only nodes with external
inputs (pinned nodes) or only nodes without external inputs (non-pinned
nodes). No cluster contains both pinned and non-pinned nodes. It follows that
either $\mathcal{C}_{k}\subset\mathcal{V}_{P}$ or
$\mathcal{C}_{k}\subset\mathcal{V}_{NP}$, $k=1,...,M$.

In order to find the equitable clusters, we apply the algorithm developed by
Belykh and Hasler belykh2011 to the extended network with Laplacian matrix
$\tilde{L}$. The extended network can always be partitioned into $M+1$
equitable clusters in which one of the clusters necessarily corresponds to the
only source node. In what follows we are going to neglect this one cluster and
focus on the remaining $M$ equitable clusters
$\mathcal{C}_{1},\mathcal{C}_{2},\ldots,\mathcal{C}_{M}$. Now we can define
the $N\times M$ equitable cluster indicator matrix $E$, where $N$ is the
number of nodes and $M$ is the number of clusters: $e_{ij}=1$ if node $i$ is
in cluster $\mathcal{C}_{j}$ and $0$ otherwise for $i=1,\ldots,N$ and
$j=1,\ldots,M$.

######  Definition 8.

Quotient Graph. Given a network with inputs, represented by the matrix pair
$(L,R)$, and the equitable clusters, the quotient network $(L_{Q},R_{Q})$ can
be defined as a network where all the nodes in each equitable cluster are
mapped to a single quotient node,

$\displaystyle L_{Q}$ $\displaystyle=(E^{T}E)^{-1/2}E^{T}LE(E^{T}E)^{-1/2}$
(30a) $\displaystyle R_{Q}$
$\displaystyle=(E^{T}E)^{-1/2}E^{T}RE(E^{T}E)^{-1/2}.$ (30b)

We introduce the orthogonal transformation matrix ${T}_{q}$ whose first $M$
rows are $T_{q}^{E}=(E^{T}E)^{-1/2}E^{T}$ and whose remaining $N-M$ rows are
orthogonal to the first ones. We call $T_{q}^{r}$ the matrix composed of the
last $N-M$ rows of the matrix ${T}_{q}$. We take the matrix $T_{q}^{r}$ to
have a particular structure, where each one of its rows is associated with one
cluster, in the sense that all of the entries in that row not corresponding to
the nodes in that cluster are equal to zero. Several such matrices have been
proposed in the literature NC ; schaub2016graph ; cho2017stable .

######  Theorem 3.

Application of the ‘quotient’ transformation ${T}_{q}$ to the pair $L,R$
yields,

$\displaystyle{T}_{q}L{T}_{q}^{-1}$ $\displaystyle=L_{Q}\oplus L_{O}$ (31a)
$\displaystyle{T}_{q}R{T}_{q}^{-1}$ $\displaystyle=R_{Q}\oplus R_{O},$ (31b)

where the square quotient matrices $L_{Q},R_{Q}$ have size $M\times M$ and the
square ‘redundant’ matrices $L_{O},R_{O}$ have size $(N-M)\times(N-M)$.

###### Proof.

For succinctness we omit the proof of (31a). We just note that this proof
follows directly from Refs. NC ; schaub2016graph ; cho2017stable .

To prove Eq. (31b), we know from Eq. (11) that the matrix $R$ is diagonal with
the first $s$ entries being one and the remaining $N-s$ entries being zero.
Therefore, the matrix $R$ has $s$ eigenvalues equal one and $N-s$ eigenvalues
equal zero. Then, we can write the matrix $R=V_{r}D_{r}V_{r}^{T}$, where the
matrix $D_{r}=D1\oplus D{0}$ such that $D_{1}$ is the $s$-dimensional identity
matrix and $D_{0}=0_{N-s}$. The matrix $V_{r}$ can be written in the form,

$V_{r}=\begin{pmatrix}({V}_{1})_{s\times s}&\mathbf{0}_{s\times(N-s)}\\\
\mathbf{0}_{(N-s)\times s}&({V}_{0})_{(N-s)\times(N-s)}\end{pmatrix}$ (32)

due to the property that the eigenvectors associated with the $1$ eigenvalue
($0$ eigenvalue) have all of their last $N-s$ (first $s$) entries equal to
zero. Now we can write,

${T}_{q}R{T}_{q}^{-1}={T}_{q}V_{r}D_{r}V_{r}^{T}{T}_{q}^{T}=({T}_{q}V_{r})D_{r}({T}_{q}V_{r})^{T}.$
(33)

It follows that

${T}_{q}R{T}_{q}^{-1}=\begin{pmatrix}T_{q}^{E_{p}}(T_{q}^{E_{p}})^{T}&T_{q}^{E_{p}}(T_{q}^{r_{p}})^{T}\\\
T_{q}^{r_{p}}(T_{q}^{E{p}})^{T}&T_{q}^{r_{p}}(T_{q}^{r_{p}})^{T}\\\
\end{pmatrix},$ (34)

where $T_{q}^{E_{p}}$ is the $M\times s$ sub-matrix that shows the restriction
of the first $m$ rows of $T_{q}$ corresponding to the clusters that contain
pinned nodes. $T_{q}^{r_{p}}$ is the $(N-M)\times s$ sub-matrix that shows the
restriction of the remaining $N-M$ rows of $T_{q}$ corresponding to the
clusters that contain non-pinned nodes. Now considering the particular
structure of $T_{q}$ and the characteristic of the equitable clusters of a
network with inputs (Remark 6), we have $T_{q}^{E_{p}}(T_{q}^{r_{p}})^{T}=0$.
Therefore, it can be concluded that ${T}_{q}R{T}_{q}^{-1}=R_{Q}\oplus R_{O}$
where $R_{Q}=T_{q}^{E_{p}}(T_{q}^{E_{p}})^{T}$ and
$R_{O}=T_{q}^{r_{p}}(T_{q}^{r_{p}})^{T}$. ∎

We conclude that application of the transformation $T_{q}$ to the pair
$(L,R)$, $T_{q}^{-1}LT_{q}$ and $T_{q}^{-1}RT_{q}$, simultaneously block
diagonalizes $(L,R)$ into an $M$-dimensional pair $(L_{Q},R_{Q})$ and an
$(N-M)$-dimensional pair $L_{O},R_{O}$.

######  Remark 7.

The results of Theorem 3 apply to any equitable cluster partition associated
with the pair $(L,R)$, not just the equitable cluster partition with the
smallest number of clusters returned by the algorithm developed by Belykh and
Hasler belykh2011 .

###  V.3 Transformation $\hat{T}$

Next we discuss the application of both transformations $T_{c}$ and $T_{q}$ to
the pair $(L,R)$. While there are different ways in which this can be done, we
found the most straightforward approach is to first apply $T_{q}$ and decouple
the pair $(L,R)$ into the pairs $(L_{Q},R_{Q})$ and $(L_{O},R_{O})$ and then
apply equivalent transformations to the pairs $(L_{Q},R_{Q})$ and
$(L_{O},R_{O})$ into their controllable and uncontrollable parts. We call
$\hat{T}$ the resulting transformation.

######  Lemma 5.

The pair $(L,R)$ is decoupled by application of the resulting transformation
$\hat{T}$ into four blocks,

$\begin{array}[]{clc}\hat{T}L\hat{T}^{-1}=\begin{pmatrix}L_{qc}&0&0&0\\\
0&L_{rc}&0&0\\\ 0&0&L_{qu}&0\\\ 0&0&0&L_{ru}\end{pmatrix}\\\ \\\
\hat{T}R\hat{T}^{-1}=\begin{pmatrix}R_{qc}&0&0&0\\\ 0&R_{rc}&0&0\\\ 0&0&0&0\\\
0&0&0&0\end{pmatrix},\end{array}$ (35)

where $(L_{qc},R_{qc})$ is the controllable quotient pair, $(L_{rc},R_{rc})$
is the controllable redundant pair, $(L_{qu},0)$ is the uncontrollable
quotient pair, and $(L_{ru},0)$ is the uncontrollable redundant pair. The size
of the pair $(L_{qc},R_{qc})$ is equal to the dimension $u$ of the
controllable subspace of the quotient pair $(L_{Q},R_{Q})$, the sum of the
sizes of the two pairs $(L_{qu},0)$ and $(L_{ru},0)$ is equal to the dimension
$N-c$ of the uncontrollable subspace of the pair $(L,R)$. The size of the pair
$(L_{rc},R_{rc})$ is equal to $c-u$.

###### Proof.

From Theorem 3 we see that application of the transformation ${T}_{q}$ to the
pair $(L,R)$ decouples the pair into a quotient pair $(L_{Q},R_{Q})$ and a
redundant pair $(L_{O},R_{O})$. Application of the equivalent transformation
to the pair $(L_{Q},R_{Q})$ will necessarily decouple the pair $(L_{Q},R_{Q})$
in its controllable and uncontrollable pairs and application of the equivalent
transformation to the pair $(L_{O},R_{O})$ will necessarily decouple the pair
$(L_{O},R_{O})$ in its controllable and uncontrollable pairs.

∎

Table 1 summarizes the application of transformation $\hat{T}$ that decouples
the pair $(L,R)$ into four block pairs, which in what follows we will refer to
as the qc pair (quotient controllable), the qu pair (quotient uncontrollable),
the rc pair (redundant controllable), and the ru pair (redundant
uncontrollable.)

Table 1: The four decoupled blocks and their dimensions after applying
transformation $\hat{T}$ to the pair ($L,R$). The first two rows show the
controllable pairs, $qc$ and $rc$, with dimensions $u$ and $c-u$,
respectively. The second two rows indicate the uncontrollable pairs, $qu$ and
$ru$. The sum of the size of the uncontrollable pairs is equal to $N-c$.

######  Theorem 4.

The transformation $\hat{T}$ leads to a finest SBD of the pair (L,R).

###### Proof.

The proof is based on the concept of common invariant subspaces
bischer2021simultaneous ; arapura2004common . As stated above, we first apply
$T_{q}$ and decouple the pair $(L,R)$ into the pairs $(L_{Q},R_{Q})$ and
$(L_{O},R_{O})$. We initially focus on the pair
$L_{Q}=T_{q}^{E}L(T_{q}^{E})^{T},R_{Q}=T_{q}^{E}R(T_{q}^{E})^{T}$ and on the
corresponding column subspaces $\mathcal{L}_{Q}$ and $\mathcal{R}_{Q}$. We
compute the matrix $T_{c}$ as described at the end of Sec. V.1 and we obtain:

$T_{c}^{-1}T_{q}^{E}L(T_{q}^{E})^{T}T_{c}=T_{c}^{-1}L_{Q}T_{c}=\begin{bmatrix}L_{qc}&0\\\
0&L_{qu}\end{bmatrix}$ (36)
$T_{c}^{-1}T_{q}^{E}R(T_{q}^{E})^{T}T_{c}=T_{c}^{-1}R_{Q}T_{c}=\begin{bmatrix}R_{qc}&0\\\
0&0\end{bmatrix}$ (37)

The column subspaces of the matrices $L_{qc}$, $R_{qc}$ and $L_{qu}$ are
$\mathcal{L}_{qc}$, $\mathcal{R}_{qc}$ and $\mathcal{L}_{q}$, respectively.
The linearly independent columns of $B_{L}=L(T_{q}^{E})^{T}$ form a basis for
the span of $L_{Q}$. Analogously, the linear independent columns of
$B_{R}=R(T_{q}^{{E}})^{T}$ form a basis for the span of $R_{Q}$. Therefore,
$\hat{B}=B_{R}^{T}B_{L}$ contains the $p$ linearly independent columns of the
product $L_{Q}R_{Q}$.

We compute the SVD $\hat{B}=\hat{U}\hat{\Sigma}\hat{V}^{T}$. The matrix
$\hat{U}_{p}$ contains the first $p$ columns of $\hat{U}$, which form a basis
of $\mathcal{R}_{Q}\cap\mathcal{L}_{Q}$ van1996matrix .

By using the basis $U_{p}$, the intersection space
$\mathcal{R}_{Q}\cap\mathcal{L}_{Q}$ can be expressed as the direct sum of
minimal-dimension subspacesgalantai2008subspaces .

The controllability matrix of the pair $(L_{Q},R_{Q})$ is
$K=[R_{Q},L_{Q}R_{Q},\dots,L_{Q}^{N-1}R_{Q}]$ and certainly contains the
columns of $\hat{B}$. Therefore, the columns of $U_{p}$ are a subset of those
of $U_{c}$ (see. Remark 5 and following text).

The matrix $T_{c}$ decomposes the vector space
$\mathcal{R}_{Q}+\mathcal{L}_{Q}$ into the following direct sum of subspaces:

$\begin{array}[]{cc}\mathcal{R}_{Q}+\mathcal{L}_{Q}=(\mathcal{L}_{qc}+\mathcal{R}_{qc})\oplus\mathcal{L}_{qu}=\\\
\\\
\left[(\mathcal{R}_{qc}\cap\mathcal{L}_{qc})\oplus(\mathcal{R}_{qc}-\mathcal{R}_{qc}\cap\mathcal{L}_{qc})\oplus(\mathcal{L}_{qc}-\mathcal{R}_{qc}\cap\mathcal{L}_{qc})\right]\oplus\mathcal{L}_{qu}\end{array}$

(38)

As stated above, $(\mathcal{R}_{qc}\cap\mathcal{L}_{qc})$ is the direct sum of
minimal-dimension subspaces because it is represented through the basis formed
by the columns of $U_{c}$. The vector space
$(\mathcal{R}_{qc}-\mathcal{R}_{qc}\cap\mathcal{L}_{qc})\oplus(\mathcal{L}_{qc}-\mathcal{R}_{qc}\cap\mathcal{L}_{qc})$
cannot be decomposed in the direct sum of smaller common invariant subspaces
because it is the direct sum of two subspaces with empty intersection.

The vector space $\mathcal{L}_{qu}$ is the direct sum of $1$-dimensional
subspaces, because it is represented through the basis formed by the columns
of $X_{N-c}$, which contain the eigenvectors of $L_{Q}$ that belong to the non
controllable subspace.

Since the transformation performed by $T_{c}$ decomposes the vector space
$\mathcal{L}_{Q}+\mathcal{R}_{Q}$ into the direct sum of minimal-dimensional
subspaces, it produces a finest SBD of the transformed matrices
$T_{c}^{-1}L_{Q}T_{c}$ e $T_{c}^{-1}R_{Q}T_{c}$.

Similar reasoning can be applied, mutatis mutandis, to the pair
$(L_{o},R_{o})$. Therefore, the transformtation $\hat{T}$ produces a finest
SBD of the pair $(L,R)$. ∎

To demonstrate the application of the transformation $\hat{T}$ to a pinning
control problem, consider the network shown in Fig. 2 (a). The network has
$M=2$ equitable clusters, one comprising only the pinned node (in yellow in
the figure) and the other one comprising all the remaining nodes (in blue in
the figure). The $M=2$-dimensional quotient network is controllable, therefore
it cannot be reduced further. By computing the transformation matrix $\hat{T}$
and following Lemma 5 and Eq. (35), we have:

$\hat{T}L\hat{T}^{-1}=\begin{pmatrix}-1.8937&2.4254&0&0&0\\\
2.4254&-3.1063&0&0&0\\\ 0&0&-1.5858&0&0\\\ 0&0&0&-3&0\\\
0&0&0&0&-4.4142\end{pmatrix}$ (39a)
$\hat{T}R\hat{T}^{-1}=\begin{pmatrix}-4.9981&0.0970&0&0&0\\\
0.0970&-0.0019&0&0&0\\\ 0&0&0&0&0\\\ 0&0&0&0&0\\\ 0&0&0&0&0\end{pmatrix}$
(39b)

By applying the transformation $\hat{T}$ to the network of Fig. 2, the pair
$(L,R)$ is decoupled into a $2$-dimensional qc block and three scalar ru
blocks. By further diagonalizing the qc block of the matrix
$\hat{T}R\hat{T}^{-1}$ and the ru block of the matrix $\hat{T}L\hat{T}^{-1}$,
we obtain the same transformation of the matrices $L$ and $R$ shown in Fig. 2
(d,e).

We also considered application of the theory to randomly generated networks.
First, we constructed matrices $L$ corresponding to $500$ Erdős-Rènyi random
networks with $N=50$ nodes and connection probability $p=0.05$. Second, to
construct the matrix $R$, we randomly selected the pinned nodes for the case
that the number $s$ was varied from $1$ to $49$. Finally, we applied both the
transformations $T$ and $\hat{T}$ to the constructed pair $(L,R)$, and we
compared the number and dimensions of the resulting blocks, which yields a
total of $24500$ numerical experiments. We repeated the process for all the
combinations of $L$ and $R$. We found that in all the cases considered, the
number and dimensions of the blocks produced by the SBD transformation $T$
matched the number and dimensions of the blocks produces by $\hat{T}$
transformation, which is in accordance with the result of Theorem 4 that
$\hat{T}$ provides a finest SBD.

We finally comment on the complexity of computing the two transformations $T$
and $\hat{T}$. The computational complexity of calculating the SBD
transformation $T$ is $O(N^{3})$ following zhang2020 . The computational
complexity of calculating the transformation $\hat{T}=T_{q}T_{c}$ is
$O(N^{3})$ kalman1963mathematical ; blondel2000survey ; cline2006computatio .
We thus conclude that the computational complexity of the two calculations is
of the same order.

##  VI Pinned node selection and multiple driven blocks

As stated in Sec. IV.2, the dimension of the driven pair $(L_{d},R_{d})$
coincides with the rank of the controllability matrix of the pair $({L},{R})$.
This indicates that both the topology of the network and the number and choice
of the pinned nodes affect the dimension of the driven pair. In general, it is
possible that the driven pair may be formed of $1\leq w\leq s$ independent
driven pairs $(L_{d}^{1},R_{d}^{1})$, $(L_{d}^{2},R_{d}^{2})$,…,
$(L_{d}^{w},R_{d}^{w})$, where $L_{d}\oplus_{i=1}^{w}L_{d}^{i}$ and
$R_{d}\oplus_{i=1}^{w}R_{d}^{i}$. If $s=1$ then there is only one driven pair.
However, the case $s>1$ in which there is more than one pinned node needs
further consideration. In what follows we show that the choice of the pinned
nodes in relation to the network topology determines $w$.

######  Definition 9.

Pinned nodes symmetry. A pinned node symmetry (PNS) is a symmetry of the pair
$(L,R)$ della2020 , i.e. a permutation matrix $\Pi$ with the following two
properties: (1) it commutes with both matrices $L$ and $R$, i.e, $L\Pi=\Pi L$,
$R\Pi=\Pi R$ and (2) it swaps two or more pinned nodes.

Because the permutation matrix $\Pi$ commutes with both $L$ and $R$, it
provides a particular solution to the problem of finding the matrix $P$
described in Sec. IV.1. Hence, the matrix $T_{\Pi}$ formed of the eigenvectors
of the matrix $\Pi$ provides an SBD transformation (though not necessarily
finest.)

######  Definition 10.

Eigenvectors of a permutation matrix garcia2015 . Each permutation matrix
$\Pi$ partitions the set of the network nodes into $f$ disjoint cycles,
$\mathcal{K}_{1},\mathcal{K}_{2},...,\mathcal{K}_{f}$ of length
$l_{1},l_{2},...,l_{f}$, respectively, $\sum_{i=1}^{f}l_{i}=N$. For instance,
a cycle of length $1$ corresponds to a node that gets mapped back to itself, a
cycle of length $2$ corresponds to two nodes that get swapped with one another
and so on. Each eigenvector of the permutation matrix is associated with just
one cycle, meaning that the entries of that eigenvector are nonzero only for
entries that correspond to nodes in the cycle. The number of eigenvectors
associated with each cycle is equal to the length of the cycle. As a result,
the matrix $T_{\Pi}$ can be written in block diagonal form,
$T_{\Pi}=\bigoplus_{i=1}^{f}T_{\Pi}^{i}$, where each square block
$T_{\Pi}^{i}$ has size $l_{i}$. For each cycle, there is one eigenvector whose
entries associated with the nodes of the cycle are all ones and the remaining
entries are all zeros. The remaining eigenvectors associated with each cycle
have the same structure and are all orthogonal to one another.

######  Lemma 6.

Choosing two or more pinned nodes that are swapped by a PNS results into
independent driven pairs for each one of these pinned nodes.

###### Proof.

The proof we provide is for the case that the PNS involves two pinned nodes.
The generalization to the case of more than two pinned nodes is
straightforward and is omitted here for conciseness. We recall that the matrix
$T_{\Pi}$ has a special structure, where each one of its rows has nonzero
entries only corresponding to the pinned nodes that are swapped by the PNS. We
can then associate a ‘quotient network’ to the PNS where all the pinned nodes
swapped by the PNS are mapped to only one node. By using Theorem 3 and Remark
7, we see that the transformation $T_{\Pi}$ decouples the pair $(L,R)$ into a
quotient block with one pinned node and an orthogonal (redundant) block with
the other pinned node. ∎

As an example, consider the $N=10$ node weighted network represented in Fig. 5
(a), with nodes color coded according to the $M=8$ equitable clusters to which
they belong. Thin (thick) edges indicate a weight of the connection equal to
one (two.) The three pinned nodes ($s=3$), number $1$, $2$, and $6$ are
indicated with a red arrow pointing at them. This is an example of a network
for which the quotient network is not fully controllable. Application of the
transformation matrix $\hat{T}$ to the pair $(L,R)$ decouples the pair $(L,R)$
into four block pairs: a $7$-dimensional qc pair (shown in black), a
$1$-dimensional rc pair (in red), a $1$-dimensional qu pair (in green) and a
$1$-dimensional ru pair (in blue).

Figure 5: (a) $N=10$-node weighted network. Nodes are colored according to the
$M=8$ equitable clusters to which they belong. Thin (thick) edges indicate a
weight of the connection equal to one (two.) Application of transformation
matrix $\hat{T}$ to (b) $L$ and (c) $R$ decouples the pair $(L,R)$ into the qc
block pair shown in black, the rc block pair shown in red, the qu block pair
shown in green, and the ru block pair shown in dark blue.

It can be seen that this network has one PNS, which corresponds to a disjoint
cycle that swaps node 1 with 2, $l=2$,

$\Pi=\begin{pmatrix}0&1&0&0&0&0&0\\\ 1&0&0&0&0&0&0\\\ 0&0&1&0&0&0&0\\\
0&0&0&1&0&0&0\\\ 0&0&0&0&1&0&0\\\ 0&0&0&0&0&1&0\\\ 0&0&0&0&0&0&1\end{pmatrix}$
(40)

The block-diagonal transformation matrix $T_{\Pi}$ whose columns are the
eigenvectors of matrix $\Pi$ is,

$T_{\Pi}=\begin{pmatrix}-1/\sqrt{2}&1/\sqrt{2}&0&0&0&0&0&0&0&0\\\
1/\sqrt{2}&1/\sqrt{2}&0&0&0&0&0&0&0&0\\\ 0&0&1&0&0&0&0&0&0&0\\\
0&0&0&1&0&0&0&0&0&0\\\ 0&0&0&0&1&0&0&0&0&0\\\ 0&0&0&0&0&1&0&0&0&0\\\
0&0&0&0&0&0&1&0&0&0\\\ 0&0&0&0&0&0&0&1&0&0\\\ 0&0&0&0&0&0&0&0&1&0\\\
0&0&0&0&0&0&0&0&0&1\end{pmatrix}$ (41)

Now by applying the transformation matrix $T_{\Pi}$ to the pair $(L,R)$, we
have:

$\begin{array}[]{clc}L_{\Pi}=\left(\leavevmode\resizebox{399.2946pt}{}{$\begin{array}[]{c | c c c c c c c c c}-3&0&0&0&0&0&0&0&0&0\\\ \hline\cr 0&-3&0&0&0&\sqrt{2}&\sqrt{2}&0&0&\sqrt{2}\\\ 0&0&-3&0&1&0&0&1&1&0\\\ 0&0&0&-3&1&0&0&1&1&0\\\ 0&0&1&1&-6&2&0&0&2&0\\\ 0&\sqrt{2}&0&0&2&-6&0&0&0&2\\\ 0&\sqrt{2}&0&0&0&0&-6&2&0&2\\\ 0&0&1&1&0&0&2&-6&2&0\\\ 0&0&1&1&2&0&0&2&-6&0\\\ 0&\sqrt{2}&0&0&0&2&2&0&0&-6\end{array}$}\right)\\\ \\\ R_{\Pi}=\left(\leavevmode\resizebox{211.3883pt}{}{$\begin{array}[]{c | c c c c c c c c c}1&0&0&0&0&0&0&0&0&0\\\ \hline\cr 0&1&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&1&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0\end{array}$}\right)\end{array}$ (42)

We see that the application of $T_{\Pi}$ to the network shown in Fig. 5,
decouples the pair $(L,R)$ into a $9$-dimensional driven quotient block and a
$1$-dimensional driven redundant block.

To conclude, Lemma 6 states that choosing two pinned nodes which are swapped
by a PNS leads to two independent driven blocks and confirms that both the
topology of the network and the number and choice of the pinned nodes affect
the number and dimensions of the driven pairs.

##  VII Application to Complex networks

In the previous sections we have applied our approach to a variety of simple
networks with a small number of nodes. Here, we extend our analysis to larger
synthetic networks and complex networks with topologies taken from the
literature. The main goal of this section is to provide evidence of what
appears to be a generic property of these networks in terms of the structure
of the blocks that arise from the decomposition. Table 2 provides for each
network information on the number of nodes $N$, the number of edges $E$, the
number of pinned nodes $s$, and $u$ the dimension of the largest block which
corresponds to the only quotient controllable block. All these networks are
undirected and unweighted. For each network in the table, $s=7$ (number chosen
without loss of generality) pinned nodes are randomly selected.

The first and second networks ($SF_{1}$ and $SF_{2}$) are synthetic scale-free
networks generated with the static model korea , with exponents of the power-
law degree distribution equal to $\alpha=2.11$ and $\alpha=2.56$,
respectively. In both cases, we first generated the network connectivity using
the static model, then removed all nodes that were not part of the network
giant component. The third dataset, DD-g147, is a biological protein
interaction network nr . The fourth dataset, ca-sandi-auths, represents a
collaboration network between scientists at Sandia National Labs nr . The last
dataset, Case60nordic is a power grid network 5491276 .

To demonstrate the use of the method to determine the stability of the
synchronous solution, we consider the Rössler oscillator (Eq. 24) at each
network node and set $\boldsymbol{G}(\boldsymbol{x}_{i})=[x_{i},y_{i},0]^{T}$
and $\boldsymbol{H}(\boldsymbol{x}_{i})=[0,y_{i},0]^{T}$. For this choice of
the functions $\boldsymbol{F},\boldsymbol{G},\boldsymbol{H}$, the target
synchronous solution is stable for $\gamma$ larger than a critical value
$\gamma^{*}_{1}$. We emphasize that the choice of the specific oscillator is
not relevant for the goal specified above, and we have obtained similar
results for other choices of the functions
$\boldsymbol{F},\boldsymbol{G},\boldsymbol{H}$. Table 2 compares the
numerically estimated values of $\gamma^{*}_{1}$ and $\gamma^{*}_{2}$, where
the latter is the critical value of $\gamma$ above which the MLE of the
largest block of the pair ($L_{T},R_{T}$) becomes negative. As can be seen,
for all these examples the MLE corresponding to the largest block determines
the stability of the target solution.

Figure 6 compares – by plotting them against the pinning control coupling
coefficient $\gamma$ – the synchronization error of the whole network (panels
a,c,e,g,i) and the MLE of the largest block of the pair ($L_{T},R_{T}$),
whereas the MLEs of the remaining blocks are all negative (panels b,d,f,h,j).

Figure 7 shows the values of $\gamma^{*}_{1}$ and $\gamma^{*}_{2}$ versus the
number of pinned nodes $s$ for the last network listed in Table 2
($60nordic$). For each value of $s$, the pinned nodes are randomly chosen. As
can be seen from Fig. 7, $\gamma_{1}^{*}$ and $\gamma_{2}^{*}$ are found to be
in close agreement for all values of $s$. Figure 7 provides evidence that,
independent of the number of pinned nodes, the stability of the target
synchronous solution is solely a function of the MLE of the largest block of
the pair $(L_{T},R_{T})$. This has two main implications: on the one hand this
is advantageous because stability of the target synchronous solution is only
determined by this one block, on the other hand it is disadvantageous because
this one block is typically quite large, indicating that our ability to reduce
the dimension of the problem is intrinsically limited for these networks.

We would like to emphasize that the purpose of this section was not to show
application of the pinning control problem to examples of practical interest,
but rather to show general properties of the blocks obtained when randomly
selecting the pin nodes of a large complex network. In this respect, we found
the remarkable property that for all the networks considered (synthetic and
‘real’) and for all the random selections of the pin nodes, there is typically
one larger quotient and controllable block, and the ability to drive the
network to the target solution is de facto determined by the stability of this
one block.

Table 2: For each network we report $N$, ${E}$, $s$, and $u$, the total number of nodes, the number of edges, the number of pinned nodes, and the dimension of the $qc$ block, respectively. $\gamma^{*}_{1}$ is the value of $\gamma$ above which the network converges to the synchronous solution, and $\gamma^{*}_{2}$ is the critical value of $\gamma$ above which the MLE of the largest block of the pair ($L_{T},R_{T}$) becomes negative. All these networks are undirected and unweighted. Label | Name | $N$ | ${E}$ | $s$ | $u$ | $\gamma^{*}_{1}$ | $\gamma^{*}_{2}$
---|---|---|---|---|---|---|---
$Net_{1}$ | $SF_{1}$ | 72 | 100 | 7 | 46 | 1.43 | 1.425
$Net_{2}$ | $SF_{2}$ | 63 | 80 | 7 | 55 | 1.43 | 1.415
$Net_{3}$ |  DD-g147nr | 99 | 317 | 7 | 94 | 2.385 | 2.368
$Net_{4}$ |  ca-sandi-authsnr | 56 | 70 | 7 | 13 | 1.27 | 1.268
$Net_{5}$ |  Case60nordiccapitanescu | 60 | 72 | 7 | 56 | 1.302 | 1.291

Figure 6: (a,c,e,g,i) Synchronization error for the pair ($L,R$) as a function
of $\gamma$. (b,d,f,h,j) MLE of the largest block of the pair ($L_{T},R_{T}$)
as a function of $\gamma$. The network labels $Net_{1}$, $Net_{2}$, etc are
consistent with Table II. Figure 7: $\gamma^{*}_{1}$ and $\gamma^{*}_{2}$
versus the number of pinned nodes for the $60nordic$ network listed in Table
2, where $\gamma^{*}_{1}$ is the value of the $\gamma$ above which the error
goes to zero, and $\gamma^{*}_{2}$ is the zero-crossing point of the MLE of
the largest block of the pair ($L_{T},R_{T}$).

##  VIII Conclusion

While techniques for the simultaneous block diagonalization (SBD) of matrices
have been previously applied to the network synchronization problem Ir:So ;
zhang2020 , here, our focus is on decoupling the synchronization stability
equations of a network with pinning control, consisting of two different types
of coupling: (i) node-to-node coupling among the network nodes and (ii) input-
to-node coupling from source node to the ‘pinned nodes’. Our main result is
that we prove that the blocks resulting from the SBD can be categorized into
four types, which we call quotient controllable, quotient uncontrollable,
redundant controllable, redundant uncontrollable. This has important
consequences as it indicates that, for this class of networks, the SBD
transformation can be replaced by application of an alternative transformation
that decouples the stability problem into a quotient and a redundant part and
each one of these two into a controllable and an uncontrollable part.
Different from previous applications of the SBD technique Ir:So ; zhang2020 ,
we provide a characterization about the dimensions of the independent sets of
equations in which the stability problem is reduced. Our analysis applied to
several complex networks from the literature shows that stability of the
target synchronous solution is always determined by the maximum Lyapunov
exponent of only one quotient controllable block.

From the standpoint of stability analysis, it is convenient to have small
blocks. One may conclude that it may be desirable to minimize the dimension of
the controllable subspace, as we have shown that the blocks corresponding to
the non-controllable subspace are all scalar, which is the best possible
outcome in terms of dimension reduction. In terms of stability, we have
investigated how the synchronizability (or pin-controllability) varies
depending on the individual blocks. One point that is left for future study is
how to choose the pinned nodes in an optimal way, i.e., so as to enforce
desirable properties of the controlled network. It is also possible to extend
our work to the case that delays affect either the node-to-node coupling
function $\boldsymbol{G}$ or the source-to-node coupling function
$\boldsymbol{H}$. The techniques for simultaneous block diagonalization of
matrices can also be applied to the case of time-varying networks
dariani2011effect ; feng2016cluster ; lin2021pinning ; shi2021synchronization
, but with the caveat that the blocks resulting from the decomposition may
vary in time with the network time evolution.

## Supplementary Material

The supplementary material includes further clarifications on the calculation
and then application of the transformation matrices $T$ and $\hat{T}$. Note 1
provides further information on the network with $N=10$ nodes shown in Figure
5. Note 2 focuses on another example of a star network with $N=4$ nodes, for
which we considered all possible choices of unique combinations of pinned
nodes. Note 3 includes the proof that that if $T^{(1)},...,T^{(M)}$ is a set
of $N\times N$ real matrices, then the intersection of their null spaces is
equal to the null space of the sum of their inner products.

## Acknowledgement

The authors are in debt to Isaac Klickstein and Galen Novello for insightful
discussions.

## Data Availability

The data that support the findings of this study are available within the
article.

##  IX Appendix

###  IX.1 Derivation of Blocks P1 and P2

Matrix-matrix products and be expressed as matrix-vector products using the
Kronecker product. Let $A\in\mathbb{R}^{n\times p}$, $X\in\mathbb{R}^{p\times
m}$ and $B\in\mathbb{R}^{n\times m}$. The matrix-matrix product $AX=B$ can be
equivalently expressed as,

$(I_{m}\otimes A)\text{vec}(X)=\text{vec}(B)$ (A.1)

Similarly, for $A\in\mathbb{R}^{p\times m}$ and $X\in\mathbb{R}^{n\times p}$
then the matrix-matrix product $XA=B$ can be equivalently expressed as,

$(A^{T}\otimes I_{n})\text{vec}(X)=\text{vec}(B)$ (A.2)

With these two identities, the four equations in Eq. (21) can be equivalently
expressed as matrix-vector product.

$((I_{s}\otimes L_{11})-(L_{11}^{T}\otimes
I_{s}))\text{vec}(P_{1})=\bm{0}_{s^{2}}$ (A.3a) $((I_{\tau}\otimes
L_{22})-(L_{22}^{T}\otimes I_{\tau}))\text{vec}(P_{2})=\bm{0}_{\tau^{2}}$
(A.3b) $(L_{12}^{T}\otimes I_{s})\text{vec}(P_{1})-(I_{\tau}\otimes
L_{12})\text{vec}(P_{2})=\bm{0}_{s\tau}$ (A.3c) $-(I_{s}\otimes
L_{21})\text{vec}(P_{1})+(L_{21}^{T}\otimes
I_{\tau})\text{vec}(P_{2})=\bm{0}_{s\tau}$ (A.3d)

Combine the first two lines of Eq. (A.3) into a composite linear system of
equations and the second two lines of Eq. (A.3) into another composite linear
system of equations to create Eq. (22).

## References

  * (1) W. Wu, W. Zhou, and T. Chen, “Cluster synchronization of linearly coupled complex networks under pinning control,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , vol. 56, no. 4, pp. 829–839, 2008.
  * (2) X. Liu and T. Chen, “Cluster synchronization in directed networks via intermittent pinning control,” _IEEE Transactions on Neural Networks_ , vol. 22, no. 7, pp. 1009–1020, 2011.
  * (3) H. Su, Z. Rong, M. Z. Chen, X. Wang, G. Chen, and H. Wang, “Decentralized adaptive pinning control for cluster synchronization of complex dynamical networks,” _IEEE Transactions on Cybernetics_ , vol. 43, no. 1, pp. 394–399, 2012.
  * (4) W. Yu, G. Chen, J. Lu, and J. Kurths, “Synchronization via pinning control on general complex networks,” _SIAM Journal on Control and Optimization_ , vol. 51, no. 2, pp. 1395–1416, 2013.
  * (5) X. Liu and T. Chen, “Synchronization of nonlinear coupled networks via aperiodically intermittent pinning control,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 26, no. 1, pp. 113–126, 2014.
  * (6) J.-L. Wang, H.-N. Wu, T. Huang, S.-Y. Ren, and J. Wu, “Pinning control for synchronization of coupled reaction-diffusion neural networks with directed topologies,” _IEEE Transactions on Systems, Man, and Cybernetics: Systems_ , vol. 46, no. 8, pp. 1109–1120, 2015.
  * (7) X. Liu and T. Chen, “Synchronization of complex networks via aperiodically intermittent pinning control,” _IEEE transactions on automatic control_ , vol. 60, no. 12, pp. 3316–3321, 2015.
  * (8) P. DeLellis, F. Garofalo, and F. L. Iudice, “The partial pinning control strategy for large complex networks,” _Automatica_ , vol. 89, pp. 111–116, 2018.
  * (9) V. N. Belykh, G. V. Osipov, V. S. Petrov, J. A. Suykens, and J. Vandewalle, “Cluster synchronization in oscillatory networks,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 18, no. 3, p. 037106, 2008\.
  * (10) V. N. Belykh, I. V. Belykh, and M. Hasler, “Hierarchy and stability of partially synchronous oscillations of diffusively coupled dynamical systems,” _Physical Review E_ , vol. 62, no. 5, p. 6332, 2000.
  * (11) G. Chen, “Pinning control of complex dynamical networks,” _IEEE Transactions on Consumer Electronics_ , 2022.
  * (12) Y. Tang, H. Gao, and J. Kurths, “Distributed robust synchronization of dynamical networks with stochastic coupling,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , vol. 61, no. 5, pp. 1508–1519, 2013\.
  * (13) G. Chen, “Pinning control and synchronization on complex dynamical networks,” _International Journal of Control, Automation and Systems_ , vol. 12, no. 2, pp. 221–230, 2014.
  * (14) W. Liu, W. Gu, W. Sheng, X. Meng, S. Xue, and M. Chen, “Pinning-based distributed cooperative control for autonomous microgrids under uncertain communication topologies,” _IEEE Transactions on Power Systems_ , vol. 31, no. 2, pp. 1320–1329, 2015.
  * (15) Y. Orouskhani, M. Jalili, and X. Yu, “Optimizing dynamical network structure for pinning control,” _Scientific reports_ , vol. 6, 2016.
  * (16) H. Chen, P. Shi, and C.-C. Lim, “Pinning impulsive synchronization for stochastic reaction–diffusion dynamical networks with delay,” _Neural Networks_ , vol. 106, pp. 281–293, 2018.
  * (17) R. Dariani, A. Buscarino, L. Fortuna, and M. Frasca, “Effect of motion on pinning control of time-varying networks,” in _Neural Nets WIRN11_. IOS Press, 2011, pp. 105–111.
  * (18) J. Feng, P. Yang, and Y. Zhao, “Cluster synchronization for nonlinearly time-varying delayed coupling complex networks with stochastic perturbation via periodically intermittent pinning control,” _Applied Mathematics and Computation_ , vol. 291, pp. 52–68, 2016.
  * (19) H. Lin and J. Wang, “Pinning control of complex networks with time-varying inner and outer coupling,” _Mathematical Biosciences and Engineering_ , vol. 18, no. 4, pp. 3435–3447, 2021.
  * (20) L. Shi, C. Zhang, and S. Zhong, “Synchronization of singular complex networks with time-varying delay via pinning control and linear feedback control,” _Chaos, Solitons & Fractals_, vol. 145, p. 110805, 2021.
  * (21) Z. An, H. Zhu, X. Li, C. Xu, Y. Xu, and X. Li, “Nonidentical linear pulse-coupled oscillators model with application to time synchronization in wireless sensor networks,” _IEEE Transactions on Industrial Electronics_ , vol. 58, no. 6, pp. 2205–2215, 2011.
  * (22) Z. Kan, J. R. Klotz, E. L. Pasiliao Jr, and W. E. Dixon, “Containment control for a social network with state-dependent connectivity,” _Automatica_ , vol. 56, pp. 86–92, 2015.
  * (23) H. L. Trentelman, K. Takaba, and N. Monshizadeh, “Robust synchronization of uncertain linear multi-agent systems,” _IEEE Transactions on Automatic Control_ , vol. 58, no. 6, pp. 1511–1523, 2013.
  * (24) G. Russo and M. D. Bernardo, “How to synchronize biological clocks,” _Journal of Computational Biology_ , vol. 16, no. 2, pp. 379–393, 2009.
  * (25) H. Su and X. Wang, _Pinning control of complex networked systems: Synchronization, consensus and flocking of networked systems via pinning_. Springer Science & Business Media, 2013.
  * (26) Y. Wu, W. Wei, G. Li, and J. Xiang, “Pinning control of uncertain complex networks to a homogeneous orbit,” _IEEE Transactions on Circuits and Systems II: Express Briefs_ , vol. 56, no. 3, pp. 235–239, 2009.
  * (27) J. G. Barajas-Ramírez, “Adaptive tracking via pinning in networks of nonidentical nodes,” _Kybernetika_ , vol. 54, no. 1, pp. 30–40, 2018.
  * (28) C. J. Vega, E. N. Sanchez, and R. Alzate, “Inverse optimal pinning control for synchronization of complex networks with nonidentical chaotic nodes,” _IFAC-PapersOnLine_ , vol. 51, no. 13, pp. 235–239, 2018.
  * (29) A. Y. Alanis, D. Ríos-Rivera, E. N. Sanchez, and O. D. Sanchez, “Learning impulsive pinning control of complex networks,” _Mathematics_ , vol. 9, no. 19, p. 2436, 2021.
  * (30) F. Uhlig, “Simultaneous block diagonalization of two real symmetric matrices,” _Linear Algebra and Its Applications_ , vol. 7, no. 4, pp. 281–289, 1973.
  * (31) T. Maehara and K. Murota, “A numerical algorithm for block-diagonal decomposition of matrix $*$-algebras with general irreducible components,” _Japan journal of industrial and applied mathematics_ , vol. 27, no. 2, pp. 263–293, 2010.
  * (32) K. Murota, Y. Kanno, M. Kojima, and S. Kojima, “A numerical algorithm for block-diagonal decomposition of matrix $*$-algebras with application to semidefinite programming,” _Japan Journal of Industrial and Applied Mathematics_ , vol. 27, no. 1, pp. 125–160, 2010.
  * (33) D. Irving and F. Sorrentino, “Synchronization of a hypernetwork of coupled dynamical systems,” _Phys. Rev. E_ , vol. 86, p. 056102, 2012.
  * (34) Y. Zhang and A. E. Motter, “Symmetry-independent stability analysis of synchronization patterns,” _SIAM Review_ , vol. 62, no. 4, pp. 817–836, 2020\.
  * (35) F. Sorrentino, M. Di Bernardo, F. Garofalo, and G. Chen, “Controllability of complex networks via pinning,” _Physical Review E_ , vol. 75, no. 4, p. 046103, 2007.
  * (36) F. Sorrentino, “Effects of the network structural properties on its controllability,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 17, no. 3, p. 033101, 2007.
  * (37) L. Pecora and T. Carroll, “Master stability functions for synchronized coupled systems,” _Phys. Rev. Lett._ , vol. 80, pp. 2109–2112, 1998.
  * (38) F. Sorrentino, “Synchronization of hypernetworks of coupled dynamical systems,” _New J. Phys._ , vol. 14, p. 033035, 2012.
  * (39) T. Maehara and K. Murota, “Algorithm for error-controlled simultaneous block-diagonalization of matrices,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 32, no. 2, pp. 605–620, 2011.
  * (40) K. Ogata, _Modern control engineering_. Prentice hall, 2010.
  * (41) T. Kailath, _Linear systems_. Prentice-Hall Englewood Cliffs, NJ, 1980, vol. 156.
  * (42) T. Maehara and K. Murota, “Algorithm for error-controlled simultaneous block-diagonalization of matrices,” _SIAM J. Matrix Anal. Appl._ , vol. 33, no. 2, pp. 605–620, 2011.
  * (43) We computed the UPO by using the continuation tool MATcont dhooge:2003 ; dhooge2006matcont .
  * (44) B. R. Hunt and E. Ott, “Optimal periodic orbits of chaotic systems,” _Physical review letters_ , vol. 76, no. 13, p. 2254, 1996.
  * (45) R. J. Sánchez-García, “Exploiting symmetry in network analysis,” _Communications Physics_ , vol. 3, no. 1, pp. 1–15, 2020.
  * (46) M. Golubitsky, I. Stewart, and D. G. Schaeffer, _Singularities and Groups in Bifurcation Theory: Volume II_. Springer Science & Business Media, 2012, vol. 69.
  * (47) M. Golubitsky and I. Stewart, _The symmetry perspective: from equilibrium to chaos in phase space and physical space_. Springer Science & Business Media, 2003, vol. 200.
  * (48) S. Martini, M. Egerstedt, and A. Bicchi, “Controllability analysis of multi-agent systems using relaxed equitable partitions,” _International Journal of Systems, Control and Communications_ , vol. 2, no. 1-3, pp. 100–121, 2010.
  * (49) I. Belykh and M. Hasler, “Mesoscale and clusters of synchrony in networks of bursting neurons,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 21, no. 1, p. 016106, 2011.
  * (50) L. M. Pecora, F. Sorrentino, A. M. Hagerstrom, T. E. Murphy, and R. Roy, “Cluster synchronization and isolated desynchronization in complex networks with symmetries,” _Nature Communications_ , vol. 5, 2014.
  * (51) M. T. Schaub, N. O’Clery, Y. N. Billeh, J.-C. Delvenne, R. Lambiotte, and M. Barahona, “Graph partitions and cluster synchronization in networks of oscillators,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 26, no. 9, p. 094821, 2016.
  * (52) Y. S. Cho, T. Nishikawa, and A. E. Motter, “Stable chimeras and independently synchronizable clusters,” _Physical review letters_ , vol. 119, no. 8, p. 084101, 2017.
  * (53) I. Bischer, C. Döring, and A. Trautner, “Simultaneous block diagonalization of matrices of finite order,” _Journal of Physics A: Mathematical and Theoretical_ , vol. 54, no. 8, p. 085203, 2021.
  * (54) D. Arapura and C. Peterson, “The common invariant subspace problem: an approach via gröbner bases,” _Linear algebra and its applications_ , vol. 384, pp. 1–7, 2004.
  * (55) C. F. Van Loan and G. Golub, _Matrix Computations_. The Johns Hopkins University Press, 1996, third Ed., Sec. 12.4.4.
  * (56) A. Galántai, “Subspaces, angles and pairs of orthogonal projections,” _Linear and Multilinear Algebra_ , vol. 56, no. 3, pp. 227–260, 2008.
  * (57) R. E. Kalman, “Mathematical description of linear dynamical systems,” _Journal of the Society for Industrial and Applied Mathematics, Series A: Control_ , vol. 1, no. 2, pp. 152–192, 1963.
  * (58) V. D. Blondel and J. N. Tsitsiklis, “A survey of computational complexity results in systems and control,” _Automatica_ , vol. 36, no. 9, pp. 1249–1274, 2000.
  * (59) A. K. Cline and I. S. Dhillon, _Computation of the Singular Value Decomposition_. CRC Press, jan 2006.
  * (60) F. Della Rossa, L. Pecora, K. Blaha, A. Shirin, I. Klickstein, and F. Sorrentino, “Symmetries and cluster synchronization in multilayer networks,” _Nature communications_ , vol. 11, no. 1, pp. 1–17, 2020.
  * (61) M. I. García Planas, M. Planas, and M. dels Dolors, “Eigenvectors of permutation matrices,” _Advances in Pure Mathematics_ , vol. 5, pp. 390–394, 2015.
  * (62) K.-I. Goh, B. Kahng, and D. Kim, “Universal behavior of load distribution in scale-free networks,” _Phys. Rev. Lett._ , vol. 87, p. 278701, 2001.
  * (63) R. A. Rossi and N. K. Ahmed, “The network data repository with interactive graph analytics and visualization,” in _AAAI_ , 2015. [Online]. Available: https://networkrepository.com
  * (64) R. D. Zimmerman, C. E. Murillo-Sánchez, and R. J. Thomas, “Matpower: Steady-state operations, planning, and analysis tools for power systems research and education,” _IEEE Transactions on Power Systems_ , vol. 26, no. 1, pp. 12–19, 2011.
  * (65) F. Capitanescu, “Suppressing ineffective control actions in optimal power flow problems,” _IET Generation, Transmission & Distribution_, vol. 14, no. 13, pp. 2520–2527, 2020.
  * (66) A. Dhooge, W. Govaerts, and Y. A. Kuznetsov, “Matcont: a matlab package for numerical bifurcation analysis of odes,” _ACM Transactions on Mathematical Software (TOMS)_ , vol. 29, no. 2, pp. 141–164, 2003.
  * (67) A. Dhooge, W. Govaerts, Y. A. Kuznetsov, W. Mestrom, A. Riet, and B. Sautois, “Matcont and cl matcont: Continuation toolboxes in matlab,” _Universiteit Gent, Belgium and Utrecht University, The Netherlands_ , 2006\.