# Cluster Synchronization of Networks via a Canonical Transformation for
Simultaneous Block Diagonalization of Matrices

Shirin Panahi  Mechanical Engineering Department, University of New Mexico,
Albuquerque, NM 87131  Isaac Klickstein  Mechanical Engineering Department,
University of New Mexico, Albuquerque, NM 87131  Francesco Sorrentino
<EMAIL_ADDRESS>Mechanical Engineering Department, University of New Mexico,
Albuquerque, NM 87131

###### Abstract

We study cluster synchronization of networks and propose a canonical
transformation for simultaneous block diagonalization of matrices that we use
to analyze stability of the cluster synchronous solution. Our approach has
several advantages as it allows us to: (1) decouple the stability problem into
subproblems of minimal dimensionality while preserving physically meaningful
information; (2) study stability of both orbital and equitable partitions of
the network nodes and (3) obtain a parametrization of the problem in a small
number of parameters. For the last point, we show how the canonical
transformation decouples the problem into blocks that preserve key physical
properties of the original system. We also apply our proposed algorithm to
analyze several real networks of interest, and we find that it runs faster
than alternative algorithms from the literature.
Keywords: Dynamical Network; Simultaneous Block Diagonalization ; Cluster
Synchronization.

The problem of cluster synchronization of networks has been studied in a
number of papers in the literature, see e.g., Pecora  _et al._ (2014);
Sorrentino  _et al._ (2016); Siddique  _et al._ (2018); Zhang and Motter
(2020), among others. This paper follows up on these previous studies and
proposes a canonical simultaneous block diagonalization routine to decouple
the stability problem into subproblems (‘blocks’) of minimal dimensionality.
Our approach has mainly two advantages: (i) each block in which the original
problem is decoupled has a clear physical interpretation and (ii) it is faster
than other algorithms proposed in the literature when applied to the analysis
of real network topologies. It also nicely reconnects with previous work
Pecora  _et al._ (2014); Sorrentino  _et al._ (2016); Siddique  _et al._
(2018) as the results in our paper are found to coincide with those in these
other papers, though the techniques used to decouple the stability problem are
different.

##  I Introduction

Cluster synchronization (CS) in networks of coupled oscillators has been the
subject of vast research efforts, see e.g., Belykh  _et al._ (2008);
Sorrentino and Ott (2007); Dahms, Lehnert, and Schöll (2012); Fu  _et al._
(2014); Kanter  _et al._ (2011); Rosin  _et al._ (2013); Williams  _et al._
(2013); Nicosia  _et al._ (2013); Schaub  _et al._ (2016). This occurs when
the network nodes are divided in clusters such that the nodes in each cluster
synchronize on the same time evolution but these time evolutions are different
for nodes in different clusters. Recent work Nicosia  _et al._ (2013); Pecora
_et al._ (2014) has elucidated the relation between the symmetries of the
network topology and the formation of clusters of synchronized dynamical units
in the network. Reference Pecora  _et al._ (2014) analyzed the formation and
stability of synchronized clusters that correspond to the orbits of the
network symmetry group. References Sorrentino  _et al._ (2016); Siddique  _et
al._ (2018); Zhang and Motter (2020) extended this study to the more general
case of equitable clusters, where the nodes in each cluster are not
necessarily symmetric, but receive the same total input from the nodes in each
one of the clusters. Cluster synchroniation in directed networks was recently
addressed in Lodi, Sorrentino, and Storace (2021).

The master stability function (MSF) approach Pecora and Carroll (1998) has
been successfully applied to characterize stability of the complete
synchronous solution for networks of coupled systems with Laplacian
connectivity. The approach is based on decoupling the stability problem into a
number of lower-dimensional problems, where each of the lower-dimensional
problems depends on an eigenvalue $p_{i}$ of the Laplacian matrix. The main
advantage of this approach is that stability of the lower dimensional system
can be parametrized in a generic parameter $p$ and one can determine the range
$R$ of the parameter $p$ over which the master stability function
$\mathcal{M}(p)<0$. Then stability can be directly assessed, _for any network_
of interest, by verifying that the relevant eigenvalues $p_{i}$ belong to the
range $R$, i.e., $p_{i}\in R$. There are, therefore, two main components of
this approach: the first one is the dimensionality reduction and the second
one is the parametrization.

The case of cluster synchronization (CS) with non-Laplacian connectivity is
more complex. In what follows, we will: (i) reduce the dimensionality of the
problem and (ii) obtain a parametrization of the lower dimensional problem in
a minimal number of parameters. The novelty of our work lies especially in
point (ii). We will see that for the case of CS we can generate a ‘canonical
transformation’ of the stability problem that corresponds to a minimal number
of parameters.

##  II PROBLEM DEFINITION: Cluster synchronization

A network of coupled dynamical systems can be described by the following set
of equations,

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

where $\boldsymbol{x}_{i}(t)$ represents the $m$-dimensional state vector of
node $i$ and $\boldsymbol{F}:R^{m}\rightarrow R^{m}$ describes the time
evolution of each individual system located at node $i$. The adjacency matrix
$A$ describes the network connectivity, i.e., $A_{ij}=A_{ji}=1$ if there is a
connection between nodes $i$ and $j$ and $A_{ij}=A_{ji}=0$ otherwise. The
function $\boldsymbol{H}:R^{m}\rightarrow R^{m}$ is the node-to-node coupling
function. We call $\mathcal{V}=\\{1,...,N\\}$ the set of the network nodes.

######  Definition 1.

Equitable cluster partition. Given the adjacency matrix $A$, representing the
network topology, one can partition the set of the network nodes $\mathcal{V}$
into subsets that we call equitable clusters,
$\mathcal{C}_{1},\mathcal{C}_{2},..,\mathcal{C}_{C}$,
$\cup_{k=1}^{C}\mathcal{C}_{k}=\mathcal{V}$,
$\mathcal{C}_{k}\cap\mathcal{C}_{\ell}=\emptyset$ for $k\neq\ell$, where

$\sum_{h\in\mathcal{C}_{\ell}}A_{ih}=\sum_{h\in\mathcal{C}_{\ell}}A_{jh},\quad\begin{aligned}
\forall i,j\in\mathcal{C}_{k}\\\
\forall\mathcal{C}_{k},\mathcal{C}_{\ell}\subset\mathcal{V}.\end{aligned}$ (2)

We call $|\mathcal{C}_{k}|=n_{k}$ the number of nodes in cluster $k=1,...,C$,
$\sum_{k=1}^{C}n_{k}=N$.

In order to find the equitable clusters, we apply the algorithm developed by
Belykh and Hasler Belykh and Hasler (2011) to the network with adjacency
matrix $A$. The algorithm returns a set of $C$ equitable clusters
$\mathcal{C}_{1},\mathcal{C}_{2},\ldots,\mathcal{C}_{C}$. Information about
each equitable cluster is contained in the $N\times N$ diagonal indicator
matrix $E_{k}=\\{E_{k_{ij}}\\}$, $k=1,..,C$, where the entry $(i,i)$ of the
matrix $E_{k}$ is equal to $1$ if node $i$ is in cluster $\mathcal{C}_{k}$ and
is equal to $0$ otherwise.

######  Definition 2.

Colored Network. The previously defined equitable clusters induce a colored
network, where each node $i$ is assigned a color $k$ if node $i$ is in cluster
$\mathcal{C}_{k}$.

Given an equitable partition of the network nodes, we can define an invariant
subspace for the set of Eqs. (1), which we call the cluster synchronization
manifold. The dynamics on this manifold is the flow-invariant cluster
synchronous time evolution Golubitsky and Stewart (2005)
$\\{\boldsymbol{s}_{1}(t),\boldsymbol{s}_{2}(t),\cdots,\boldsymbol{s}_{C}(t)\\}$,
where $\boldsymbol{s}_{1}(t)$ is the synchronous solution for all nodes in
cluster $\mathcal{C}_{1}$, $\boldsymbol{s}_{2}(t)$ is the synchronous solution
for nodes in cluster $\mathcal{C}_{2}$, and so on.

The $C\times C$ quotient matrix $Q$ is defined such that for each pair of
equitable clusters $\mathcal{C}_{k}$ and $\mathcal{C}_{l}$ we have,

$Q_{kl}=\sum_{j\in\mathcal{C}_{l}}A_{ij}\quad i\in\mathcal{C}_{k}.$ (3)

The quotient matrix describes a network (the ‘quotient network’), for which
all the nodes in each equitable cluster collapse to a single quotient node. By
assuming the system of equations (1) evolves on the cluster synchronization
manifold, and averaging over all the nodes in each cluster, we can derive the
equations for the time evolution of the quotient network,

$\dot{\boldsymbol{s}}_{k}(t)=\boldsymbol{F}(\boldsymbol{s}_{k}(t))+\sum_{l=1}^{C}Q_{kl}\boldsymbol{H}(\boldsymbol{s}_{v}(t)),\quad
k,l=1,2,\cdots,C,$ (4)

where the $m$-dimensional vector $\boldsymbol{s}_{k}(t)$ represents the state
of the quotient network node $k=1,...,C$.

######  Definition 3.

Equitable clusters encoding matrix. With knowledge of the clusters, we can
construct the $N\times C$ equitable cluster encoding matrix $O$, such that
$O_{ij}=1$ if node $i$ is in cluster $\mathcal{C}_{j}$ and $0$ otherwise for
$i=1,\ldots,N$ and $j=1,\ldots,C$.

######  Remark 1.

Given a network described by the adjacency matrix $A$ and the equitable
clusters encoding matrix $O$ (see definition 3), the $C\times C$-dimensional
quotient matrix $Q$ can be computed as follows Schaub  _et al._ (2016),

$Q=(O^{T}O)^{-1}O^{T}AO.$ (5)

To investigate the stability of the cluster synchronous solution, we consider
a small perturbation
$\delta\boldsymbol{x}_{i}=(\boldsymbol{x}_{i}-\boldsymbol{s}_{k})$,
$i\in\mathcal{C}_{k}$. By linearizjng Eq. (1) about Eq. (4) we obtain the
vectorial equation,

$\delta\dot{\boldsymbol{x}}(t)=\left[\sum_{c=1}^{C}E_{c}\otimes
D\boldsymbol{F}(\boldsymbol{s}_{c}(t))+\sum_{c=1}^{C}AE_{c}\otimes
D\boldsymbol{H}(\boldsymbol{s}_{c}(t))\right]\delta\boldsymbol{x}(t),$ (6)

in the $mN$\- dimensional vector
$\delta{\boldsymbol{x}}(t)=[\delta{\boldsymbol{x}}_{1}^{T}(t),\delta{\boldsymbol{x}}_{2}^{T}(t),...,\delta{\boldsymbol{x}}_{N}^{T}(t)]^{T}$.
We are interested in the possibility that the stability problem for the
$mN$-dimensional system Eq. (6) can be decoupled into a set of lower-
dimensional equations. To this end, we will show how to obtain a
transformation matrix $T$ that block-diagonalizes the matrices $AE_{c}$ while
leaving the matrices $E_{c}$ unchanged.

##  III Dimensionality Reduction of the Cluster Synchronization Stability
Problem

As mentioned in the previous section, our goal is to reduce the stablility
problem into a set of independent lower-dimensional equations instead of
dealing with the high-dimensional problem, Eq. (6). Ref. Pecora  _et al._
(2014) has used a transformation based on the irreducible representations of
the symmetry group to block-diagonalize the set of Eqs. (6) for the case of a
network with symmetries. A limitation of this approach is that it is only
applicable to the case of ‘orbital clusters’ (corresponding to the orbits of
the symmetry group, see e.g. Tinkham (2003)) and not to the more general case
of ‘equitable clusters’ Kudose (2009). Stability of the cluster synchronous
solution in the case of equitable clusters was first discussed in Refs.
Sorrentino  _et al._ (2016); Siddique  _et al._ (2018). An alternative method
based on the approach of simultaneous block diagonalization of matrices (SBD)
Maehara and Murota (2010a, b, 2011), firstly applied to the problem of
complete synchronization of networks in Irving and Sorrentino (2012), was used
for cluster synchronization in Zhang and Motter (2020). Different from Zhang
and Motter (2020), in this paper we propose an SBD transformation that will
lead to a ‘canonical SBD transformation’ of the stability problem. Our method
applies to both the cases of equitable and orbital clusters. We will discuss
the benefits of this canonical transformation in the rest of this paper.

The problem of simultaneous block diagonalization can be formalized as
follows: given a set of $N\times N$ symmetric matrices $A^{(1)},...,A^{(M)}$
find an $N\times N$ orthogonal matrix $T$ such that the matrices
$T^{-1}A^{(k)}T$ have a common block-diagonal structure for $k=1,...,M$. It
should be noted that such a block-diagonal structure is not unique in at least
two senses: first, the blocks may be permuted, resulting in block diagonal
decompositions that are isomorphic; second, the matrices corresponding to
certain blocks may be further refined into smaller blocks, resulting in a
structures that are fundamentally different. A block diagonal structure with
smaller blocks is considered to be finer and the finest SBD is beneficial as
it provides the simplest elements in the decoupling of systems as described
above.

In the rest of this paper, we will always focus on finding a finest SBD. We
write

$T=\mathcal{SBD}(A^{(1)},A^{(2)},...,A^{(M)})$ (7)

to indicate that the following transformation yields

$T^{-1}A^{(k)}T=B^{(k)},\quad k=1,...,M,$ (8)

where all the matrices $B^{(k)}$, $k=1,...,M$ share the same finest block
diagonal form,

$B^{(k)}=\bigoplus_{j}B_{j}^{(k)},$ (9)

with the blocks $B_{j}^{(k)}$ all having the same sizes $\beta_{k}$ for
$k=1,...,M$ and not being further reducible by a simultaneous transformation.

The method for calculating a finest SBD transformation proposed in Maehara and
Murota (2011) requires two main steps: (i) Finding a matrix $P$ that commutes
with a set of matrices, and (ii) calculating the transformation matrix $T$
formed of the eigenvectors of the matrix $P$.

###  III.1 Matrix P

In order to study stability of the cluster synchronization solution, we need
to simultaneously block diagonalize the set of $C+1$ matrices
$\\{A,E_{1},E_{2},\cdots,E_{C}\\}$ where the matrix $A$ is the adjacency
matrix of the network and each $E_{i}$ is the binary and diagonal cluster
indicator matrix Zhang and Motter (2020). Once we find
$T=\mathcal{SBD}(A,E_{1},E_{2},...,E_{C})$, the matrices $A$ and $E_{i}$ are
transformed as follows,

$T^{-1}AT=B=\oplus_{k=1}^{r}\hat{B}^{k},$ (10a)
$T^{-1}E_{i}T=J_{i}=\oplus_{k=1}^{r}\hat{J}^{k}_{i}\quad i=1,\cdots,C,$ (10b)

where $B$ is the transformed matrix $A$, $J_{i}$ is the transformed matrix
$E_{i}$, and the blocks of the $C+1$-tuple
$(\hat{B}^{k},\hat{J}^{k}_{1},....,\hat{J}^{k}_{C})$ have the same dimensions
$\beta_{k}$, $k=1,...,r$, $\sum_{k=1}^{r}\beta_{k}=N$. In what follows, we
will refer to $(\hat{B}^{k},\hat{J}^{k}_{1},....,\hat{J}^{k}_{C})$ as a block
tuple.

Without loss of generality, we assume that the network nodes are ordered so
that the first $n_{1}$ nodes are the ones in cluster $\mathcal{C}_{1}$,
followed by the $n_{2}$ nodes in cluster $\mathcal{C}_{2}$, and so on, and the
last $n_{C}$ nodes are the ones in cluster $\mathcal{C}_{C}$. Then the
matrices $E_{k}$ are in the following form,

$E_{1}=\begin{pmatrix}I_{n_{1}}&0\\\ 0&{0}_{N-n_{1}}\end{pmatrix}\quad
E_{2}=\begin{pmatrix}{0}_{n_{1}}&0&0\\\ 0&I_{n_{2}}&0\\\
0&0&{0}_{N-(n_{1}+n_{2})}\end{pmatrix}\quad\cdots\quad
E_{C}=\begin{pmatrix}{0}_{N-n_{C}}&0\\\ 0&I_{n_{C}}\end{pmatrix}$ (11)

where $I_{n_{i}}$ is the identity matrix of size $n_{i}$, and $0_{n_{i}}$ is
the zero matrix of size $n_{i}$.

######  Lemma 1.

Any matrix $P$ that commutes with the set of matrices $E_{i}$ in Eq. (11) has
the following block-diagonal structure,

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

where the block $P_{1}$ has dimension $n_{1}$ and the block $P_{2}$ has
dimension $n_{2}$ and so on.

###### Proof.

The product $PE_{i}$ must have columns of zeros in the same place as $E_{i}$,
so if $P$ is to commute with $E_{i}$ then $E_{i}P$ must also have these
columns as zero. Because $E_{i}$ has an identity block in its non-zero
columns, this forces $P$ to have zeros in the intersection rows where $E_{i}$
has its identity block and the columns where $E_{i}$ is zero. Placing all such
zeros as necessary for each of the $E_{i}$ gives the structure of $P$ as
claimed. ∎

Following Lemma 1, we can see that the transformation matrix $T$, which has
the eigenvectors of the matrix $P$ for its columns, must also have the same
block-diagonal structure as the matrix $P$,

$T=\begin{pmatrix}T_{1}&0&\cdots&0\\\ 0&T_{2}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ 0&0&\cdots&T_{C}\end{pmatrix}$ (13)

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

######  Remark 2.

A trivial consequence of Eq. (13) is that as the matrix $T$ is orthogonal,
each one of the block matrices $T_{k}$ is also orthogonal,
$T_{k}^{T}=T_{k}^{-1}$.

######  Remark 3.

The particular block-diagonal structure of the matrix $T$, where each block
$T_{k}$ corresponds to cluster $\mathcal{C}_{k}$ has important consequences.
Each indipendent block obtained from application of the SBD transformation
will only produce linear combinations of the states of nodes from the same
cluster (same color). This implies that to each block corresponds a colored
subnetwork (corresponding to the block), see Definition 2. It also follows
that each block describes stability of either one cluster or a set of
intertwined clusters Pecora  _et al._ (2014); Siddique  _et al._ (2018).

######  Remark 4.

An additional benefit of our matrix $T$ is that each one of its rows is
associated to one and only one of the clusters, in the sense that all the
entries that do not correspond to the nodes of that cluster are zero. This is
important as it indicates the particular way in which the cluster is broken if
the Lyapunov exponent of the corresponding block of the matrix $B$ is positive
Siddique  _et al._ (2018), i.e., entries that are the same (different)
correspond to nodes that remain (do not remain) synchronized after the
breaking.

######  Lemma 2.

Application of the block diagonal transformation matrix $T$ to each cluster
indicator matrix $E_{i}$ ensures that $J_{i}=E_{i}$ for $i=1,2,\cdots,C$. In
other words, each $E_{i}$ gets mapped back to itself using the transformation
matrix $T$ of Eq. (13).

###### Proof.

By considering the matrices $\\{E_{1},E_{2},\cdots,E_{C}\\}$ of Eq. (11) and
the block diagonal transformation matrix $T$ of Eq. (13) and by using the fact
that a block-diagonal matrix can be inverted block by block, we have
$T^{-1}E_{i}T=E_{i}$ for $i=1,2,\cdots,C$. For instance, for $E_{1}$ we have:

$T^{-1}E_{1}T=\begin{pmatrix}T^{-1}_{1}&0&\cdots&0\\\ 0&T^{-1}_{2}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&T^{-1}_{C}\end{pmatrix}\begin{pmatrix}I_{n1}&0&\cdots&0\\\
0&0&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&0\end{pmatrix}\begin{pmatrix}T_{1}&0&\cdots&0\\\
0&T_{2}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&T_{C}\end{pmatrix}=\begin{pmatrix}I_{c1}&0&\cdots&0\\\
0&0&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\ 0&0&\cdots&0\end{pmatrix}$ (14)

∎

######  Lemma 3.

Each block $T_{i}$ of the matrix $T$ has a column with entries that are all
the same.

###### Proof.

According to Lemma 1 and Eq. (13), we know that both the matrices $P$ and $T$
are block-diagonal in $C$ blocks, where each block corresponds to one cluster.
Also, according to Ref. Sánchez-García (2020) the adjacency matrix $A$
corresponding to a network with $C$ equitable clusters, has $C$ ‘non-
redundant’ eigenvectors $\mathbf{v}_{k}$, each one associated with a cluster
$\mathcal{C}_{k}$, $k=1,...,C$ whose entries ${\mathbf{v}_{k}}_{i}$ are all
the same for $i\in\mathcal{C}_{k}$ and zero otherwise. From the facts that the
matrix $P$ commutes with $A$ and has block-diagonal structure (Eq.
$\eqref{Pblock}$) it follows that the matrix $P$ will also have eigenvectors
$\mathbf{v}_{k}$, $k=1,...,C$. Thus each block $T_{k}$ associated with cluster
$\mathcal{C}_{k}$ has a column with entries that are all the same and equal to
${n_{k}}^{-1/2}$. ∎

######  Remark 5.

The transformation matrix $T$ has two properties in common with the
transformation introduced in Refs. Pecora  _et al._ (2014); Della Rossa  _et
al._ (2020): one is the block diagonal structure (13) and the other one is the
property (Lemma 3) that each block has a column with entries that are all the
same. Thus the two transformations are similar, with the following
differences: (1) the transformation matrix $T$ is not based on calculation of
the network symmetries and can be successfully applied to the case of
equitable clusters and (2) the transformation matrix $T$ is faster to compute.

######  Lemma 4.

One of the block tuples resulting from the SBD transformation of the set
$\\{A,E_{1},E_{2},...,E_{C}\\}$ is $C$-dimensional and corresponds to the
quotient network dynamics.

###### Proof.

As the $T$ matrix is orthogonal $T^{-1}=T^{T}$. From the matrix $T^{T}$ it is
possible to extract $C$ rows, where each row $i=1,...,C$ includes the column
of the block $T_{i}$ with entries that are all the same and equal to
$n_{i}^{-1/2}$. We call $\Delta$ the matrix obtained by stacking together
these $C$ rows. We note that the so-constructed matrix
$\Delta=(O^{T}O)^{-1/2}O^{T}$. Then the $C\times C$-dimensional block,

$\hat{B}^{1}=\Delta A\Delta^{T}=(O^{T}O)^{1/2}Q(O^{T}O)^{-1/2},$ (15)

from which we see that the two matrices $\hat{B}^{1}$ and $Q$ are similar. A
similar proof can be found in Klickstein, Pecora, and Sorrentino (2019).

∎

######  Remark 6.

Lemma 4 is useful as it allows us to identify one block-tuple that is
associated with dynamics parallel to cluster synchronization manifold (the
dynamics of the quotient network.) We use the label $k=1$ to indicate this
block-tuple $(\hat{B}^{1},\hat{J}^{1}_{1},....,\hat{J}^{1}_{C})$. We use the
labels $k=2,...,C$ to indicate the block-tuples associated with dynamics
transverse to the synchronization manifold. The transverse block-tuples
determine stability of the cluster synchronous solution Pecora  _et al._
(2014).

######  Lemma 5.

The sub-matrices $P_{i}\quad i=1,2,\cdots,C$ can be found by computing the
null subspace of the $\sum_{i}n_{i}^{2}\times\sum_{i}n_{i}^{2}$-dimensional
matrix $S^{T}S$ defined below.

###### Proof.

By rewriting the adjacency matrix $A$ as:

$\begin{pmatrix}A_{11}&A_{12}&\cdots&A_{1C}\\\ A_{21}&A_{22}&\cdots&A_{2C}\\\
\vdots&\vdots&\ddots&\vdots\\\ A_{C1}&A_{C2}&\cdots&A_{CC}\\\ \end{pmatrix}$
(16)

the commutation equation $P{A}={A}P$ becomes

$\begin{pmatrix}P_{1}&0&\cdots&0\\\ 0&P_{2}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&P_{C}\end{pmatrix}\begin{pmatrix}A_{11}&A_{12}&\cdots&A_{1C}\\\
A_{21}&A_{22}&\cdots&A_{2C}\\\ \vdots&\vdots&\ddots&\vdots\\\
A_{C1}&A_{C2}&\cdots&A_{CC}\\\
\end{pmatrix}=\begin{pmatrix}A_{11}&A_{12}&\cdots&A_{1C}\\\
A_{21}&A_{22}&\cdots&A_{2C}\\\ \vdots&\vdots&\ddots&\vdots\\\
A_{C1}&A_{C2}&\cdots&A_{CC}\\\ \end{pmatrix}\begin{pmatrix}P_{1}&0&\cdots&0\\\
0&P_{2}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&P_{C}\end{pmatrix}$ (17a)
$\begin{pmatrix}P_{1}A_{11}&P_{1}A_{12}&\cdots&P_{1}A_{1C}\\\
P_{2}A_{21}&P_{2}A_{22}&\cdots&P_{2}A_{2C}\\\ \vdots&\vdots&\ddots&\vdots\\\
P_{C}A_{C1}&P_{C}A_{C2}&\cdots&P_{C}A_{CC}\\\
\end{pmatrix}=\begin{pmatrix}A_{11}P_{1}&A_{12}P_{2}&\cdots&A_{1C}P_{C}\\\
A_{21}P_{1}&A_{22}P_{2}&\cdots&A_{2C}P_{C}\\\ \vdots&\vdots&\ddots&\vdots\\\
A_{C1}P_{1}&A_{C2}P_{2}&\cdots&A_{CC}P_{C}\\\ \end{pmatrix}.$ (17b)

Equation (17b) corresponds to the following $C^{2}$ equations which should all
be simultaneously solved:

$P_{i}A_{ij}-A_{ij}P_{j}=0_{n_{i},n_{j}}\quad i,j=1,2,\cdots,C$ (18)

where $n_{i}$ is the dimension of the sub-matrix $P_{i}$.

Define the function $\text{vec}:\mathbb{R}^{n\times m}\mapsto\mathbb{R}^{nm}$
that maps a matrix to a vector by stacking the matrix’s columns. We obtain the
set of $C^{2}$ equations,

$(A_{ij}^{T}\otimes I_{n_{pi}})vec(P_{i})-(I_{n_{pi}}\otimes
A_{ij})vec(P_{j})=0_{n_{i},n_{j}}$ (19)

Accordingly, Eq. (19) can be expressed as two linear systems of equations in
the matrices,

$S_{1}=\begin{pmatrix}\bar{A}_{1}&0&\cdots&0\\\ 0&\bar{A}_{2}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ 0&0&\cdots&\bar{A}_{C}\end{pmatrix}$ (20a)
$S_{2}=\begin{pmatrix}A_{12}^{T}\otimes I_{n_{1}}&-I_{n_{2}}\otimes
A_{12}&0&0&\cdots&0&0\\\ -I_{n_{1}}\otimes A_{21}&A_{21}^{T}\otimes
I_{n_{2}}&0&0&\cdots&0&0\\\ A_{13}^{T}\otimes I_{n_{1}}&0&-I_{n_{3}}\otimes
A_{13}&0&\cdots&0&0\\\ -I_{n_{1}}\otimes A_{31}&0&A_{31}^{T}\otimes
I_{n_{3}}&0&\cdots&0&0\\\ \vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
A_{1C}^{T}\otimes I_{n_{1}}&0&0&0&\cdots&0&-I_{n_{C}}\otimes A_{1C}\\\
-I_{n_{1}}\otimes A_{C1}&0&0&0&\cdots&0&A_{C1}^{T}\otimes I_{n_{C}}\\\
0&A_{23}^{T}\otimes I_{n_{2}}&-I_{n_{3}}\otimes A_{23}&0&\cdots&0&0\\\
0&-I_{n_{2}}\otimes A_{32}&A_{32}^{T}\otimes I_{n_{3}}&0&\cdots&0&0\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\ 0&A_{2C}^{T}\otimes
I_{n_{2}}&0&0&\cdots&0&-I_{n_{C}}\otimes A_{2C}\\\ 0&-I_{n_{2}}\otimes
A_{C2}&0&0&\cdots&0&A_{C2}^{T}\otimes I_{n_{C}}\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&0&0&\cdots&A_{C-1,C}^{T}\otimes I_{n_{(C-1)}}&-I_{n_{C}}\otimes
A_{C-1,C}\\\ 0&0&0&0&\cdots&-I_{n_{(C-1)}}\otimes
A_{C,C-1}&A_{C,C-1}^{T}\otimes I_{n_{C}}\end{pmatrix}$ (20b)

where $\bar{A}_{i}=(A_{ii}^{T}\otimes I_{n_{i}})-(I_{n_{i}}\otimes A_{ii})$.
By stacking Eq. (20a) and (20b) together as
$S=\begin{bmatrix}S_{1}^{T}&S_{2}^{T}\end{bmatrix}^{T}$, we search for a
vector
$\bm{p}=\begin{bmatrix}\text{vec}(P_{1})^{T}&\text{vec}(P_{2})^{T}&\cdots&\text{vec}(P_{C})^{T}\end{bmatrix}^{T}$
that lies in $\mathcal{N}(S)$, the nullspace of $S$. The matrix
$S\in\mathbb{R}^{N_{r}\times N_{c}}$ has dimensions
$N_{r}=\sum_{i=1}^{C}\sum_{j=1}^{C}n_{i}n_{j}$ is the number of rows, and
$N_{c}=\sum_{i=1}^{C}n_{i}^{2}$ is the number of columns. By using the
property of nullspaces that $\mathcal{N}(S)=\mathcal{N}(S^{T}S)$, we instead
look for a vector $\bm{p}$ that lies in $\mathcal{N}(S^{T}S)$,

$S^{T}S\bm{p}=\bm{0}_{N_{c}},$ (21)

where $S^{T}S\in\mathbb{R}^{N_{c}\times N_{c}}$. ∎

######  Remark 7.

It is worth mentioning that the original method Maehara and Murota (2011)
requires finding a vector in the nullspace of a $N^{2}\times
N^{2}$-dimensional matrix while the method derived here requires finding a
vector in the nullspace of a
$\sum_{i}n_{i}^{2}\times\sum_{i}n_{i}^{2}$-dimensional matrix where
$\sum_{i}n_{i}^{2}$ is not greater than $N^{2}$.

$n^{2}=(\sum_{i}n_{i})^{2}\geq\sum_{i}n_{i}^{2}$ (22)

Equality only occurs if $C=1$ while the minimum of $\sum_{i}n_{i}^{2}$ is
achieved when $C=N$ (so $n_{i}=1$).

By applying $T$ to Eq. (6), we obtain,

$\begin{array}[]{lcl}\dot{\boldsymbol{\eta}}(t)=\left[(T^{-1}\sum_{k=1}^{C}E_{k}T)\otimes
D\boldsymbol{F}(\boldsymbol{s}_{k}(t))+(T^{-1}\sum_{k=1}^{C}AE_{k}T)\otimes
D\boldsymbol{H}(\boldsymbol{s}_{k}(t))\right]\boldsymbol{\eta}(t),\end{array}$
(23)

where $\boldsymbol{\eta}(t)=(T^{-1}\otimes I_{m})\delta\boldsymbol{x}(t)$. Eq.
(23) can be rewritten,

$\dot{\boldsymbol{\eta}}(t)=\left[\sum_{k=1}^{C}J_{k}\otimes
D\boldsymbol{F}(\boldsymbol{s}_{k}(t))+\sum_{k=1}^{C}BJ_{k}\otimes
D\boldsymbol{H}(\boldsymbol{s}_{k}(t))\right]{\boldsymbol{\eta}}(t),$ (24)

which can be decoupled into $r$ independent equations of smaller dimensions,
in the blocks of the block diagonal matrices $J_{k}$ and $BJ_{k}$.

###  III.2 Importance of the canonical transformation

We call an SBD transformation of the set of $C+1$ matrices
$\\{A,E_{1},E_{2},...,E_{C}\\}$ ‘canonical’ if the transformation maps the
matrices $E_{1},E_{2},...,E_{C}$ back to themselves. There are at least two
main strengths of a canonical SBD transformation: (i) the parametrization
associated with the SBD transformation, and (ii) the interpretation of the
matrices resulting from this transformation. We briefly comment on both
strenghts below. Consider that the blocks obtained by the SBD transformation
of the set of $C+1$ matrices $\\{A,E_{1},E_{2},...,E_{C}\\}$ have sizes
$\beta_{k}$, $k=1,...,r$, $\sum_{k=1}^{r}\beta_{k}=N$. Then the number of
nonzero entries that parametrize the matrices after the transformation is
equal to $p_{1}=(C+1)\sum_{k=1}^{r}\beta_{k}(\beta_{k}+1)/2$. In the case of a
canonical SBD transformation the number of nonzero entries that parametrize
the matrices after the transformation is equal to
$p_{2}=\sum_{k=1}^{r}\beta_{k}(\beta_{k}+1)/2$, $p_{2}<p_{1}$, taking into
account that the matrices $J_{1},J_{2},...,J_{C}$ are known ‘a priori’. Also
in terms of interpretation, the diagonal matrices $J_{1},J_{2},...,J_{C}$ with
either zeros or ones on the main diagonal indicate that each network node
after the transformation belongs to one and only one cluster, i.e., node $i$
belongs (does not belong) to cluster $k$ if the entry $(i,i)$ of matrix
$J_{k}$ is one (zero.) Therefore, one can still say that the resulting ‘nodes’
have a color (corresponding to the cluster they belong to) while node-node
interactions are limited to the transformed matrix $B$, see also remarks 3 and
4. This is different from the case of a non-canonical transformation where the
resulting ‘nodes’ do not necessarily have a color (i.e., they do not
correspond to any of the clusters of the original network) and there are in
general $C+1$ ‘layers of connectivity’ between these nodes.

The SBD transformation described in Sec. III.1 is canonical while the SBD
transformation proposed in Zhang and Motter (2020) is not. For example,
consider the $N=4$-dimensional network with $C=2$ clusters, shown in Fig. 1
(a). The nodes are colored according to the cluster to which they belong.

Figure 1: (a) An $N=4$ node network with $C=2$ clusters. Nodes are color-coded
according to the clusters to which they belong. (b) The adjacency matrix of
the network ($A$) and the two cluster indicator matrices ($(E_{1},E_{2})$) are
graphically shown, in which a nonzero entry is indicated as a black dot. (c)
and (e) show the quotient and transverse sub-networks after application of the
SBD transformation $T$ and $\tilde{T}$, respectively. Subnetwork nodes are
colored based on the cluster with which they are associated. Black nodes
indicate that they are not associated with only one cluster. In (e)
connections with different colors represent coupling through the $\bar{B}$
subnetwork and the $\bar{J}_{i}$ subnetworks. (d) The set of matrices
$\\{B,J_{1},J_{2}\\}$ and (f) the set of matrices
$\\{\bar{B},\bar{J}_{1},\bar{J}_{2}\\}$ obtained by using the transformation
$T$ and $\tilde{T}$, respectively. The background color of each matrix entry
indicates the cluster about which the dynamics is linearized. A gray
background color indicates that the dynamics is linearized about multiple
clusters.

Figure 1 shows the results after applying a canonical (c, d) and a non-
canonical Zhang and Motter (2020) (e, f) SBD transformation. It can be seen
that both transformations reduce the $4m$-dimensional stability problem of Eq.
(6) into a $2m$-dimensional equation corresponding to the quotient dynamics
and a $2m$-dimensional equation corresponding to the transverse dynamics. Only
the latter is responsible for stability of the cluster synchronous solution.
Equation (6) for the network of Fig. 1 is:

$\begin{array}[]{cc}\delta\dot{{\bf
X}}=\begin{pmatrix}\delta\dot{\boldsymbol{x}_{1}}\\\
\delta\dot{\boldsymbol{x}_{2}}\\\ \delta\dot{\boldsymbol{x}_{3}}\\\
\delta\dot{\boldsymbol{x}_{4}}\\\
\end{pmatrix}=\begin{pmatrix}{\color[rgb]{1,0.49,0}D\boldsymbol{F}(\boldsymbol{s}_{1})}&0&0&0\\\
0&{\color[rgb]{1,0.49,0}D\boldsymbol{F}(\boldsymbol{s}_{1})}&0&0\\\
0&0&{\color[rgb]{0,0.6,0}D\boldsymbol{F}(\boldsymbol{s}_{2})}&0\\\
0&0&0&{\color[rgb]{0,0.6,0}D\boldsymbol{F}(\boldsymbol{s}_{2})}\\\
\end{pmatrix}\delta{\bf X}\\\
+\begin{pmatrix}0&0&0&{\color[rgb]{0,0.6,0}D\boldsymbol{H}(\boldsymbol{s}_{2})}\\\
0&0&{\color[rgb]{0,0.6,0}D\boldsymbol{H}(\boldsymbol{s}_{2})}&0\\\
0&{\color[rgb]{1,0.49,0}DH(\boldsymbol{s}_{1})}&0&{\color[rgb]{0,0.6,0}D\boldsymbol{H}(\boldsymbol{s}_{2})}\\\
{\color[rgb]{1,0.49,0}D\boldsymbol{H}(\boldsymbol{s}_{1})}&0&{\color[rgb]{0,0.6,0}D\boldsymbol{H}(\boldsymbol{s}_{2})}&0\\\
\end{pmatrix}\delta{\bf X}.\end{array}$ (25)

Using the canonical SBD transformation matrix $T$, Eq. (25) becomes,

$\begin{array}[]{cc}\dot{{\bf Y}}=T^{-1}\delta\dot{{\bf
X}}=\frac{\sqrt{2}}{2}\begin{pmatrix}\delta\dot{\boldsymbol{x}}_{1}+\delta\dot{\boldsymbol{x}}_{2}\\\
-\delta\dot{\boldsymbol{x}}_{3}-\delta\dot{\boldsymbol{x}}_{4}\\\
\delta\dot{\boldsymbol{x}}_{1}-\delta\dot{\boldsymbol{x}}_{2}\\\
\delta\dot{\boldsymbol{x}}_{3}-\delta\dot{\boldsymbol{x}}_{4}\\\
\end{pmatrix}=\begin{pmatrix}{\color[rgb]{1,0.49,0}D\boldsymbol{F}(\boldsymbol{s}_{1})}&0&0&0\\\
0&{\color[rgb]{0,0.6,0}D\boldsymbol{F}(\boldsymbol{s}_{2})}&0&0\\\
0&0&{\color[rgb]{1,0.49,0}D\boldsymbol{F}(\boldsymbol{s}_{1})}&0\\\
0&0&0&{\color[rgb]{0,0.6,0}D\boldsymbol{F}(\boldsymbol{s}_{2})}\\\
\end{pmatrix}{\bf Y}\\\
+\begin{pmatrix}0&-{\color[rgb]{0,0.6,0}D\boldsymbol{H}(\boldsymbol{s}_{2})}&0&0\\\
-{\color[rgb]{1,0.49,0}D\boldsymbol{H}(\boldsymbol{s}_{1})}&{\color[rgb]{0,0.6,0}D\boldsymbol{H}(\boldsymbol{s}_{2})}&0&0\\\
0&0&0&{\color[rgb]{0,0.6,0}D\boldsymbol{H}(\boldsymbol{s}_{2})}\\\
0&0&{\color[rgb]{1,0.49,0}D\boldsymbol{H}(\boldsymbol{s}_{1})}&-{\color[rgb]{0,0.6,0}D\boldsymbol{H}(\boldsymbol{s}_{2})}\\\
\end{pmatrix}{\bf Y},\end{array}$ (26)

where the vectors ${\bf Y}=[{\bf Y}_{1}^{T},{\bf Y}_{2}^{T},{\bf
Y}_{3}^{T},{\bf Y}_{4}^{T}]^{T}$ in Eq. (26) is the transformed vector
$\delta{\bf X}$ using canonical ($T$) transformation. We see that the dynamics
of each transformed variable ${\bf Y}_{i}$ is described by either one of the
two Jacobians $D\boldsymbol{F}(\boldsymbol{s}_{1})$ and
$D\boldsymbol{F}(\boldsymbol{s}_{2})$. This is due to the fact that the
canonical transformation $T$ only produces linear combinations of nodes
belonging to the same cluster (see Remark 3.) For example, the yellow nodes in
Fig. 1 (a) get mapped back to other yellow nodes in Fig. 1 (c), and so on.
These transformed nodes are then coupled through the connections of only one
network corresponding to the resulting matrix $B$ (in two blocks.)

Using the non-canonical SBD transformation matrix $\tilde{T}$, Eq. (25)
becomes,

${\begin{array}[]{ccc}\dot{{\bf Z}}=\tilde{T}^{-1}\delta\dot{{\bf
X}}=\begin{pmatrix}-0.65\delta\dot{\boldsymbol{x}}_{1}+0.65\delta\dot{\boldsymbol{x}}_{2}-0.26\delta\dot{\boldsymbol{x}}_{3}+0.26\delta\dot{\boldsymbol{x}}_{4}\\\
-0.26\delta\dot{\boldsymbol{x}}_{1}+0.26\delta\dot{\boldsymbol{x}}_{2}+0.65\delta\dot{\boldsymbol{x}}_{3}-0.65\delta\dot{\boldsymbol{x}}_{4}\\\
0.69\delta\dot{\boldsymbol{x}}_{1}+0.69\delta\dot{\boldsymbol{x}}_{2}-0.16\delta\dot{\boldsymbol{x}}_{3}-0.16\delta\delta\dot{\boldsymbol{x}}_{4}\\\
0.16\delta\dot{\boldsymbol{x}}_{1}+0.16\delta\dot{\boldsymbol{x}}_{2}+0.69\delta\dot{\boldsymbol{x}}_{3}+0.69\delta\dot{\boldsymbol{x}}_{4}\\\
\end{pmatrix}\\\
=\begin{pmatrix}0.14{\color[rgb]{1,0.49,0}D\boldsymbol{F}(\boldsymbol{s}_{1})}+0.86{\color[rgb]{0,0.6,0}D\boldsymbol{F}(\boldsymbol{s}_{2})}&0.35{\color[rgb]{1,0.49,0}D\boldsymbol{F}(\boldsymbol{s}_{1})}-0.35{\color[rgb]{0,0.6,0}D\boldsymbol{F}(\boldsymbol{s}_{2})}&0&0\\\
0.35{\color[rgb]{1,0.49,0}D\boldsymbol{F}(\boldsymbol{s}_{1})}-0.35{\color[rgb]{0,0.6,0}D\boldsymbol{F}(\boldsymbol{s}_{2})}&0.86{\color[rgb]{1,0.49,0}D\boldsymbol{F}(\boldsymbol{s}_{1})}+0.14{\color[rgb]{0,0.6,0}D\boldsymbol{F}(\boldsymbol{s}_{2})}&0&0\\\
0&0&0.05{\color[rgb]{1,0.49,0}D\boldsymbol{F}(\boldsymbol{s}_{1})}+0.95{\color[rgb]{0,0.6,0}D\boldsymbol{F}(\boldsymbol{s}_{2})}&0.22{\color[rgb]{1,0.49,0}D\boldsymbol{F}(\boldsymbol{s}_{1})}-0.22{\color[rgb]{0,0.6,0}D\boldsymbol{F}(\boldsymbol{s}_{2})}\\\
0&0&0.22{\color[rgb]{1,0.49,0}D\boldsymbol{F}(\boldsymbol{s}_{1})}-0.22{\color[rgb]{0,0.6,0}D\boldsymbol{F}(\boldsymbol{s}_{2})}&0.95{\color[rgb]{1,0.49,0}D\boldsymbol{F}(\boldsymbol{s}_{1})}+0.05{\color[rgb]{0,0.6,0}D\boldsymbol{F}(\boldsymbol{s}_{2})}\\\
\end{pmatrix}{\bf Z}\\\
+\begin{pmatrix}0.35{\color[rgb]{1,0.49,0}D\boldsymbol{H}(\boldsymbol{s}_{1})}-0.51{\color[rgb]{0,0.6,0}D\boldsymbol{H}(\boldsymbol{s}_{2})}&0.86{\color[rgb]{1,0.49,0}D\boldsymbol{H}(\boldsymbol{s}_{1})}+0.20{\color[rgb]{0,0.6,0}D\boldsymbol{H}(\boldsymbol{s}_{2})}&0&0\\\
-0.14{\color[rgb]{1,0.49,0}D\boldsymbol{H}(\boldsymbol{s}_{1})}+1.2{\color[rgb]{0,0.6,0}D\boldsymbol{H}(\boldsymbol{s}_{2})}&-0.35{\color[rgb]{1,0.49,0}D\boldsymbol{H}(\boldsymbol{s}_{1})}-0.49{\color[rgb]{0,0.6,0}D\boldsymbol{H}(\boldsymbol{s}_{2})}&0&0\\\
0&0&0.22{\color[rgb]{1,0.49,0}D\boldsymbol{H}(\boldsymbol{s}_{1})}+1.17{\color[rgb]{0,0.6,0}D\boldsymbol{H}(\boldsymbol{s}_{2})}&0.95{\color[rgb]{1,0.49,0}D\boldsymbol{H}(\boldsymbol{s}_{1})}-0.27{\color[rgb]{0,0.6,0}D\boldsymbol{H}(\boldsymbol{s}_{2})}\\\
0&0&-0.05{\color[rgb]{1,0.49,0}D\boldsymbol{H}(\boldsymbol{s}_{1})}+0.72{\color[rgb]{0,0.6,0}D\boldsymbol{H}(\boldsymbol{s}_{2})}&-0.22{\color[rgb]{1,0.49,0}D\boldsymbol{H}(\boldsymbol{s}_{1})}-0.17{\color[rgb]{0,0.6,0}D\boldsymbol{H}(\boldsymbol{s}_{2})}\\\
\end{pmatrix}{\bf Z},\end{array}}$ (27)

where the vector ${\bf Z}=[{\bf Z}_{1}^{T},{\bf Z}_{2}^{T},{\bf
Z}_{3}^{T},{\bf Z}_{4}^{T}]^{T}$ in Eq. (27) is the transformed vector
$\delta{\bf X}$ using the non-canonical ($\tilde{T}$) SBD transformation. We
see that the transformation matrix $\tilde{T}$ produces the same number of
blocks with the same dimensions as the transformation matrix $T$. However,
from Eq. (27) and Fig. 1 we see that application of the non-canonical SBD
transformation does not preserve the color of the nodes in the original
network. Also, we see that the individual ${\bf Z}_{i}$ are coupled with one
another through both the $D\boldsymbol{F}$ and the $D\boldsymbol{H}$ terms.
This is shown in Fig. 1(e) where the resulting nodes are colored black to
indicate that they do not correspond to any individual color of the original
network. Connections of different colors are used to indicate the different
types of coupling.

As it is stated before, the stability analysis of the cluster synchronization
problem (Eq. (6)) is based on two main steps: the first one is the
dimensionality reduction and the second one is the parametrization. Here,
concerning the second step, we focus on the importance of a canonical SBD
transformation in order to obtain a parametrization of the lower dimensional
problem in a minimal number of parameters. For example, consider the network
from Ref. Pecora  _et al._ (2014) shown in Fig. 2 (a). We set
$q_{1}=A_{1,8}=A_{8,1}$ and $q_{2}=A_{5,10}=A_{10,5}$ and parametrize the
problem in terms of $q_{1}$ and $q_{2}$. The two links between nodes $1$ and
$8$ and nodes $5$ and $10$ are highlighted in yellow. Fig. 2 (b) and (c) show
the block-diagonal matrix $B$ and $\bar{B}$ obtained from application of the
canonical ($T$) and non-canonical ($\tilde{T}$) transformation, respectively.
In both cases, we indicate with $q_{1}$ ($q_{2}$) the block entries that are
affected by variations of the parameters $q_{1}$ ($q_{2}$.) As we can see,
changing $q_{1}$ and $q_{2}$ only affects a few entries of the block diagonal
matrix in the case of the canonical transformation, but affects many more
entries in the case of a non-canonical transformation. Looking at the
transverse $2$-dimensional block in (b) we see that $q_{1}$ only affects one
entry of this block and $q_{2}$ only affects one other entry of that block. On
the other hand, each entry of the $2$-dimensional block in (c) is affected by
both $q_{1}$ and $q_{2}$.

Figure 2: (a) $N=11$-dimensional network from Ref. Pecora  _et al._ (2014)
with two weighted edges. The two edges between nodes $1$ and $8$ and nodes $5$
and $10$ are highlighted in yellow. We set $q_{1}=A_{1,8}=A_{8,1}$ and
$q_{2}=A_{5,10}=A_{10,5}$ and parametrize the problem in terms of $q_{1}$ and
$q_{2}$. (b) shows the matrix $B$ obtained from application of a canonical SBD
transformation and (c) shows the matrix $\bar{B}$ obtained from application of
a non-canonical SBD transformation. In both cases, we label with $q_{1}$
($q_{2}$) the block entries that are affected by variations of the parameters
$q_{1}$ ($q_{2}$.)

##  IV Orbital and Equitable Partitions

We are now going to apply our approach to the case of a network for which the
equitable and orbital clusters do not coincide. This network is from Ref.
Kudose (2009). As for the other examples in this paper, the number and
dimensions of the blocks will be the same as those obtained from the method in
Ref. Zhang and Motter (2020) but the structure of the blocks will be
different. We will comment on the benefits of using a canonical
transformation.

We demonstrate the benefits of the canonical transformation with the
$N=8$-dimensional network from Ref. Kudose (2009) shown in Fig. 3 and 4.
Figure 3 is for the case of the network equitable partition with $C=2$
clusters and Fig. 4 is for the case of the network orbital partition with
$C=3$ clusters. We compare the results after application of the canonical SBD
transformation $T$ and the non-canonical SBD transformation $\tilde{T}$.
Figure 3 shows that in both cases the SBD leads to a total of $5$ blocks, one
$2$-dimensional block corresponding to the quotient network, one
$3$-dimensional block and three $1$-dimensional blocks, all corresponding to
the transverse dynamics.

An interesting case is that of the $3$-dimensional block, which is shown in
panels (c) for the case of the canonical transformation and (e) for the case
of the non-canonical transformation. In panel (c) the dynamics of each node is
linearized about the dynamics of one and only one of the quotient network
nodes (and so it retains the color of that node), while in panel (e) each node
is associated to a linear combination of the dynamics of the quotient network
nodes. Also the nodes in (e) are coupled through different types of
connections, consistent with the non-diagonal structure of the 3-dimensional
blocks in (f).

Another example of a network for which the equitable and orbital clusters do
not coincide from Ref. Siddique  _et al._ (2018) is presented in the
Supplementary Information. For this other example, for both the case of the
orbital and equitable partition, we obtain the same decomposition in blocks
already presented in Ref. Siddique  _et al._ (2018).

Figure 3: (a) An $N=8$-dimensional network with $C=2$ equitable clusters.
Nodes are color-coded according to the clusters to which they belong. (b) The
adjacency matrix of the network ($A$) and the two cluster indicator matrices
($(E_{1},E_{2})$) are graphically shown, in which a nonzero entry is indicated
as a black dot. (c) and (e) show the quotient and transverse sub-networks
after application of the SBD transformation $T$ and $\tilde{T}$, respectively.
Subnetwork nodes are colored based on the cluster to which they are
associated. Black nodes indicate that they are not associated with only one
cluster. In (e) connections with different colors represent coupling through
the $\bar{B}$ subnetwork and the $\bar{J}_{i}$ subnetworks. (d) The set of
matrices $\\{B,J_{1},J_{2}\\}$ and (f) the set of matrices
$\\{\bar{B},\bar{J}_{1},\bar{J}_{2}\\}$ obtained by using the transformation
$T$ and $\tilde{T}$, respectively. The background color of each matrix entry
indicates the cluster about which the dynamics is linearized. A gray
background color indicates that the dynamics is linearized about multiple
clusters. Figure 4: (a) An $N=8$-dimensional network with $C=3$ orbital
clusters. Nodes are color-coded according to the clusters to which they
belong. (b) The adjacency matrix of the network ($A$) and the three cluster
indicator matrices ($(E_{1},E_{2},E_{3})$) are graphically shown, in which a
nonzero entry is indicated as a black dot. (c) and (e) show the quotient and
transverse sub-networks after application of the SBD transformation $T$ and
$\tilde{T}$, respectively. Subnetwork nodes are colored based on the cluster
to which they are associated. Black nodes indicate that they are not
associated with only one cluster. In (e) connections with different colors
represent coupling through the $\bar{B}$ subnetwork and the $\bar{J}_{i}$
subnetworks. (d) The set of matrices $\\{B,J_{1},J_{2},J_{3}\\}$ and (f) the
set of matrices $\\{\bar{B},\bar{J}_{1},\bar{J}_{2},\bar{J}_{3}\\}$ obtained
by using the transformation $T$ and $\tilde{T}$, respectively. The background
color of each matrix entry indicates the cluster about which the dynamics is
linearized. A gray background color indicates that the dynamics is linearized
about multiple clusters.

##  V Real Networks Analysis

It is important to consider application of the algorithm developed in this
paper to real network topologies. In Table 1 we apply both transformation $T$
and $\tilde{T}$ to block diagonalize the adjacency matrices of several real
networks from the literature. For each network dataset, we include information
on the number of nodes $N$, the number of edges $E$, the number of nontrivial
clusters $N_{ntc}$ (clusters with more than one node Pecora  _et al._ (2014)),
the size of the largest equitable cluster $\max(|{n_{c}}|)$ and the runtime in
seconds for calculating the transformation matrices $T$ and $\tilde{T}$
(method from Zhang and Motter (2020)). All these networks are connected,
undirected and unweighted (only the giant component was considered in the case
of networks that are not connected.) For all the networks we have analyzed,
the two transformation matrices $T$ and $\tilde{T}$ produced equivalent SBD
decompositions, with the same number of blocks and of the same sizes. However,
the computation time of our algorithm was considerably (up to ten times)
lower. Our algorithm is also faster than those in Pecora  _et al._ (2014);
Siddique  _et al._ (2018) which are based on computing the irreducible
representations of the network symmetry group.

Table 1: Real networks analysis. $N$ is the number of nodes, $E$ the number of
edges, $N_{ntc}$ is the number of nontrivial clusters, $\max(|{n_{c}}|)$ is
the size of the largest equitable cluster. We include the runtime in seconds
for calculation of the transformation matrix $T$ and of the transformation
matrix $\tilde{T}$ using the code from Zhang and Motter (2020).
Name |  $N$ |  $E$ |  $N_{ntc}$ |  $\max(|{n_{c}}|)$ |  Runtime $\tilde{T}$Zhang and Motter (2020) |  Runtime ${T}$
---|---|---|---|---|---|---
ca-netscience: Scientist Collaboration Network Rossi and Ahmed (2015) | 379 | 914 | 70 | 6 | 55.7192 | 9.0219
Chilean Power Grid NetworkKim  _et al._ (2018); Bhatta, Hayat, and Sorrentino (2021) | 218 | 527 | 29 | 7 | 3.8168 | 1.5133
Power Grid Network of Western Germany Matke, Medjroubi, and Kleinhans (2016) | 491 | 665 | 43 | 5 | 103.2976 | 11.5019
Metabolic NetworkRossi and Ahmed (2015) | 453 | 2025 | 28 | 4 | 39.5829 | 10.0304
Us AirlineBatagelj and Mrvar (2006) | 332 | 2126 | 31 | 12 | 17.8030 | 5.7232
can-634Rossi and Ahmed (2015) | 429 | 1312 | 20 | 3 | 44.7595 | 13.2443
Erdos971Rossi and Ahmed (2015) | 634 | 3297 | 4 | 2 | 107.1852 | 33.2083

##  VI Conclusions

In this paper we have studied cluster synchronization of networks and proposed
a canonical transformation for simultaneous block diagonalization of matrices
that we use to study stability analysis of the cluster synchronous (CS)
solution. Our approach presents several advantages as it allows us to: (1)
decouple the stability problem into blocks of minimal dimensionality, while
preserving physically meaningful information; (2) study stability of the CS
solution for both the cases of orbital and equitable partitions of the network
nodes and (3) obtain a parametrization of the problem in a minimal number of
parameters.

Our algorithm is faster than others Sorrentino  _et al._ (2016); Siddique  _et
al._ (2018); Zhang and Motter (2020) when applied to real network topologies
and leads to a decomposition of the stability problem into a number of
subproblems that preserve key physical properties of the original system (such
as the ‘color’ of the nodes). The main advantage of our canonical
transformation is that it allows a parametrization of the stability problem in
a small number of parameters. In particular, we show how the stability
analysis for different networks can be studied in terms of these parameters
and how changing the coupling strengths of a subset of connections is
reflected in ‘localized’ variations of these parameters.

## Acknowledgement

The authors thank Galen Novello for insightful conversations on the subject of
$*$-algebra.

## Author Declarations

The authors have no conflicts to disclose.

## Data Availability

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

## Code Availability

The MATLAB code to compute the simultaneous block diagonalizations for the
cluster synchronization examples shown in this paper can be accessed from its
Github repo.

## References

  * Batagelj and Mrvar (2006) Batagelj, V. and Mrvar, A., “Pajek datasets (2006),” (2006).
  * Belykh and Hasler (2011) Belykh, I. and Hasler, M., “Mesoscale and clusters of synchrony in networks of bursting neurons,” Chaos: An Interdisciplinary Journal of Nonlinear Science 21, 016106 (2011).
  * Belykh  _et al._ (2008) Belykh, V. N., Osipov, G. V., Petrov, V. S., Suykens, J. A., and Vandewalle, J., “Cluster synchronization in oscillatory networks,” Chaos: An Interdisciplinary Journal of Nonlinear Science 18, 037106 (2008).
  * Bhatta, Hayat, and Sorrentino (2021) Bhatta, K., Hayat, M. M., and Sorrentino, F., “Modal decomposition of the linear swing equation in networks with symmetries,” IEEE Transactions on Network Science and Engineering (2021).
  * Dahms, Lehnert, and Schöll (2012) Dahms, T., Lehnert, J., and Schöll, E., “Cluster and group synchronization in delay-coupled networks,” Physical Review E 86, 016202 (2012).
  * Della Rossa  _et al._ (2020) Della Rossa, F., Pecora, L., Blaha, K., Shirin, A., Klickstein, I., and Sorrentino, F., “Symmetries and cluster synchronization in multilayer networks,” Nature communications 11, 1–17 (2020).
  * Fu  _et al._ (2014) Fu, C., Lin, W., Huang, L., and Wang, X., “Synchronization transition in networked chaotic oscillators: The viewpoint from partial synchronization,” Physical Review E 89, 052908 (2014).
  * Golubitsky and Stewart (2005) Golubitsky, M. and Stewart, I., “Synchrony versus symmetry in coupled cells,” in  _EQUADIFF 2003_ (World Scientific, 2005) pp. 13–24.
  * Irving and Sorrentino (2012) Irving, D. and Sorrentino, F., “Synchronization of dynamical hypernetworks: Dimensionality reduction through simultaneous block-diagonalization of matrices,” Physical Review E 86, 056102 (2012).
  * Kanter  _et al._ (2011) Kanter, I., Kopelowitz, E., Vardi, R., Zigzag, M., Kinzel, W., Abeles, M., and Cohen, D., “Nonlocal mechanism for cluster synchronization in neural circuits,” EPL (Europhysics Letters) 93, 66001 (2011).
  * Kim  _et al._ (2018) Kim, H., Olave-Rojas, D., Álvarez-Miranda, E., and Son, S.-W., “In-depth data on the network structure and hourly activity of the central chilean power grid,” Scientific data 5, 1–10 (2018).
  * Klickstein, Pecora, and Sorrentino (2019) Klickstein, I., Pecora, L., and Sorrentino, F., “Symmetry induced group consensus,” Chaos: An Interdisciplinary Journal of Nonlinear Science 29, 073101 (2019).
  * Kudose (2009) Kudose, S., “Equitable partitions and orbit partitions,” Acta Mathematica Sinica , 1–9 (2009).
  * Lodi, Sorrentino, and Storace (2021) Lodi, M., Sorrentino, F., and Storace, M., “One-way dependent clusters and stability of cluster synchronization in directed networks,” Nature Communications 12 (2021).
  * Maehara and Murota (2010a) Maehara, T. and Murota, K., “Error-controlling algorithm for simultaneous block-diagonalization and its application to independent component analysis,” JSIAM Letters 2, 131–134 (2010a).
  * Maehara and Murota (2010b) Maehara, T. and Murota, K., “A numerical algorithm for block-diagonal decomposition of matrix $\ast$-algebras with general irreducible components,” Japan journal of industrial and applied mathematics 27, 263–293 (2010b).
  * Maehara and Murota (2011) Maehara, T. and Murota, K., “Algorithm for error-controlled simultaneous block-diagonalization of matrices,” SIAM Journal on Matrix Analysis and Applications 32, 605–620 (2011).
  * Matke, Medjroubi, and Kleinhans (2016) Matke, C., Medjroubi, W., and Kleinhans, D., “SciGRID - An Open Source Reference Model for the European Transmission Network (v0.2),”  (2016).
  * Nicosia  _et al._ (2013) Nicosia, V., Valencia, M., Chavez, M., Díaz-Guilera, A., and Latora, V., “Remote synchronization reveals network symmetries and functional modules,” Physical review letters 110, 174102 (2013).
  * Pecora and Carroll (1998) Pecora, L. and Carroll, T., “Master stability functions for synchronized coupled systems,” Phys. Rev. Lett. 80, 2109–2112 (1998).
  * Pecora  _et al._ (2014) Pecora, L. M., Sorrentino, F., Hagerstrom, A. M., Murphy, T. E., and Roy, R., “Cluster synchronization and isolated desynchronization in complex networks with symmetries,” Nature communications 5, 1–8 (2014).
  * Rosin  _et al._ (2013) Rosin, D. P., Rontani, D., Gauthier, D. J., and Schöll, E., “Control of synchronization patterns in neural-like boolean networks,” Physical review letters 110, 104102 (2013).
  * Rossi and Ahmed (2015) Rossi, R. A. and Ahmed, N. K., “The network data repository with interactive graph analytics and visualization,” in  _Proceedings of the Twenty-Ninth AAAI Conference on Artificial Intelligence_ (2015).
  * Sánchez-García (2020) Sánchez-García, R. J., “Exploiting symmetry in network analysis,” Communications Physics 3, 1–15 (2020).
  * Schaub  _et al._ (2016) Schaub, M. T., O’Clery, N., Billeh, Y. N., Delvenne, J.-C., Lambiotte, R., and Barahona, M., “Graph partitions and cluster synchronization in networks of oscillators,” Chaos: An Interdisciplinary Journal of Nonlinear Science 26, 094821 (2016).
  * Siddique  _et al._ (2018) Siddique, A. B., Pecora, L., Hart, J. D., and Sorrentino, F., “Symmetry-and input-cluster synchronization in networks,” Physical Review E 97, 042217 (2018).
  * Sorrentino and Ott (2007) Sorrentino, F. and Ott, E., “Network synchronization of groups,” Phys. Rev. E 76, 056114 (2007).
  * Sorrentino  _et al._ (2016) Sorrentino, F., Pecora, L. M., Hagerstrom, A. M., Murphy, T. E., and Roy, R., “Complete characterization of stability of cluster synchronization in complex dynamical networks,” Science Advances 2 (2016).
  * Tinkham (2003) Tinkham, M.,  _Group theory and quantum mechanics_ (Courier Corporation, 2003).
  * Williams  _et al._ (2013) Williams, C. R., Murphy, T. E., Roy, R., Sorrentino, F., Dahms, T., and Schöll, E., “Experimental observations of group synchrony in a system of chaotic optoelectronic oscillators,” Physical review letters 110, 064104 (2013).
  * Zhang and Motter (2020) Zhang, Y. and Motter, A. E., “Symmetry-independent stability analysis of synchronization patterns,” SIAM Review 62, 817–836 (2020).