# Failure of the simultaneous block diagonalization technique applied to
complete and cluster synchronization of random networks

Shirin Panahi, Nelson Amaya, Isaac Klickstein, Galen Novello, Francesco
Sorrentino
University of New Mexico, Albuquerque, US 80131

###### Abstract

We discuss here the application of the algebraic techniques for the
simultaneous block diagonalization (SBD) of matrices to the study of the
stability of both complete and cluster synchronization in random (generic)
networks. For both problems, we see that the extent of the dimensionality
reduction achievable is the same as that produced by application of a trivial
transformation.

Keywords: Network; Simultaneous Block Diagonalization ; Complete
Synchronization; Cluster Synchronization.

##  I Introduction

The algebraic technique for simultaneous block diagonalization (SBD) of
matrices has been originally proposed by Maehara, Murota and collaborators in
a number of seminal papers Maehara and Murota (2010a, b); Murota  _et al._
(2010); Maehara and Murota (2011). The first paper where this technique was
applied to network synchronization is Irving and Sorrentino (2012), which
focused on complete synchronization of networks with nodes connected through
two or more coupling functions. More recently this technique has been applied
to cluster synchronization of networks Zhang and Motter (2020). The use of
this technique allows to reduce the stability problem in a number of blocks of
smallest dimension. The original Master Stability Function (MSF) derivation of
Pecora and Carroll (1998b) decouples the stability problem for any
$N$-dimensional matrix corresponding to an undirected network into $N$
independent blocks, where each block coincides with one of the matrix
eigenvalues. However, it is unclear what the extent of the dimensionality
reduction obtained from applicaton of the SBD technique can be. Here we
characterize the extent of this dimensionality reduction when the SBD approach
is applied to generic networks, where by a generic network we mean a ‘typical’
network that is produced by a random process such as the Erdős-Rényi network
generation algorithm Erdős and Rényi (1960) or the configuration model Molloy
and Reed (1995). We show that application of the SBD reduction to these random
(generic) networks does not lead to a beneficial reduction of the stability
problem, either in the case of complete synchronization or cluster
synchronization. With this, we do not mean that the technique is not useful.
However, it points out that its usefulness is limited to the non-generic case,
for which the reduction can sometimes yet be very significant Irving and
Sorrentino (2012); Zhang and Motter (2020).

##  II Simultaneous Block Diagonalization of Matrices

The problem of simultaneous block diagonalization can be formalized as
follows: given a set of $N\times N$ matrices $A^{(1)},...,A^{(M)}$ find an
$N\times N$ orthogonal matrix $P$ such that the matrices $P^{T}A^{(k)}P$ 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 (FSBD) is beneficial in that it
provides the simplest elements in the decoupling of systems as described
above.

There are two different but closely related theoretical frameworks with which
we can address our problem of finding a block-diagonal decomposition for a
finite set of given $N\times N$ real matrices. The one is group representation
theory Miller (1973); Serre (1977) which relies on group symmetries and
ensures a degree of universality in a SBD. The second is matrix $*$-algebras
Wedderburn (1934) which are not only necessary to answer the fundamental
theoretical question of the existence of such a finest block-diagonal
decomposition but also useful in its computation. Indeed, existence can be
justified through the structure theorem of $*$-algebras (Kojima, Kojima, and
Hara, 1997, Theorem 5.4) and this structure has also been utilized to
formulate algorithms for computing the SBD of $A^{(1)},...,A^{(M)}$. In
particular, our approach appeals to this structure, but it should be noted
that both frameworks have been utilized in the literature Bai  _et al._
(2009); De Klerk and Sotirov (2010); Gatermann and Parrilo (2004); Riener  _et
al._ (2013); Kanno  _et al._ (2001).

In what follows we will write

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

to indicate that the following transformation yields

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

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

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

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

##  III Procedure to Determine $P$

Here we describe the procedure to compute the FSBD for a set of $M$ symmetric
matrices denoted $A^{(k)}$, $k=1,2,\ldots,M$, previously published Maehara and
Murota (2011). First, we find a matrix $U$ that simultaneously commutes with
each of matrices $A^{(k)}$ Maehara and Murota (2011), that is,
$[A^{(k)},U]=A^{(k)}U-UA^{(k)}=O_{n}$, $k=1,\ldots,M$, where $O_{n}$ is the
$n$-by-$n$ matrix of all zeroes.

Define the vectorizing function $\text{vec}:\mathbb{R}^{n\times
m}\mapsto\mathbb{R}^{nm}$ to take as input an $n$-by-$m$ matrix and return a
vector by stacking each of the matrix’s columns on top of each into a vector
of length $nm$. For two matrices $A\in\mathbb{R}^{n\times m}$ and
$B\in\mathbb{R}^{p\times\ell}$, let $A\otimes B\in\mathbb{R}^{np\times m\ell}$
denote the Kronecker product. As the commutator equation is linear in $U$, it
can alternatively be expressed as a matrix-vector product.

$\text{vec}(A^{(k)}U-UA^{(k)})=\text{vec}(A^{(k)}U-\text{vec}(UA^{(k)})=\text{vec}(O_{n}),\quad
k=1,2,\ldots,M$ (4)

The vectorizing function applied to a matrix product can be expressed as a
matrix-vector product, where for $A\in\mathbb{R}^{n\times m}$ and
$B\in\mathbb{R}^{m\times p}$, the product $\text{vec}(AB)=(I_{p}\otimes
A)\text{vec}(B)=(B^{T}\otimes I_{n})\text{vec}(A)$ (see Proposition 7.1.9 in
Bernstein (2009)). Apply these identities to Eq. (4), and define
$P^{(k)}\text{vec}(U)=(I_{n}\otimes A^{(k)}-(A^{(k)})^{T}\otimes
I_{n})\text{vec}(U)=\bm{0}_{N^{2}}$ where $\bm{0}_{N^{2}}$ is the vector of
all zeros of length $N^{2}$. To find $U$, we look for a vector in the
intersection of the nullspaces of $P^{(k)}$ for $k=1,2,\ldots,M$, that is, a
vector $\text{vec}(U)\in\bigcap_{k=1}^{M}\mathcal{N}(P^{(k)})$. This can be
accomplished in two steps by first noting that for a matrix
$A\in\mathbb{R}^{n\times m}$ $\mathcal{N}(A)=\mathcal{N}(A^{T}A)$ (see Theorem
2.4.3 in Bernstein (2009)) and second, for a set of $M$ positive semi-definite
matrices $B^{(j)}$, $j=1,2,\ldots,M$,
$\mathcal{N}\left(\sum_{j=1}^{M}B^{(j)}\right)=\bigcap_{j=1}^{M}\mathcal{N}(B^{(j)})$
(see Fact 8.7.3 in Bernstein (2009)). As the matrix $P^{(k)}$ may not be
positive semi-definite, the vectorized commutator operation is pre-multiplied
by $(P^{(k)})^{T}$, so that the matrix $(P^{(k)})^{T}P^{(k)}$ is symmetric and
positive semi-definite. Create the matrix
$S=\sum_{k=1}^{M}(P^{(k)})^{T}P^{(k)}$ so that if a vector $\text{vec}(U)$ is
in the nullspace of $S$, it lies in the intersection of the nullspaces of
$P^{(k)}$, and thus it also commutes with all $A^{(k)}$, $k=1,\ldots,M$.
To determine the nullspace of the matrix $S\in\mathbb{R}^{NN\times NN}$, which
by construction is positive semi-definite, we find the eigenvectors of $S$
corresponding to eigenvalues equal to zero. While $S$ is large and dense, its
special structure makes finding a few extremal eigenvalues and eigenvectors
feasible even when $N$ is large by using the Lanczos method Saad (2003) which
only requires, from the user, a function to compute matrix-vector products.
Note that while a matrix-vector product, $S\bm{u}$, requires $N^{4}$
operations (remember $S\in\mathbb{R}^{NN\times NN}$), it can equivalently be
computed using nested commutation operations required $4N^{3}$ operations. To
see this reduction, we can break the matrix-vector product into individual
contributions from each $P^{(k)}$.

$S\bm{u}=\sum_{k=1}^{M}P^{(k)^{T}}P^{(k)}\bm{u}=\sum_{k=1}^{M}P^{(k)^{T}}\text{vec}\left([A^{(k)},\text{mat}(U)]\right)=\sum_{k=1}^{M}\text{vec}\left(\left[A^{(k)^{T}},\left[A^{(k)},U\right]\right]\right)$
(5)

Evaluating a commutator requires two $N$-by-$N$ matrix products which each
requires $2N^{3}$ operations. In total, for $M$ twice nested commutators, the
total work required is $4MN^{3}$ operations, which for $M\ll N$, is a
significant reduction as compared to constructing $S$ explicitly. This can be
demonstrated with the following steps to compute
$\text{vec}(Y)=S\text{vec}(U)$.

  1. 1.

Initialize $Y=O_{N}$ to be the $N$-by-$N$ matrix of zeroes

  2. 2.

For $k=1,2,\ldots,M$

    1. (a)

$\hat{U}=A^{(k)}U-UA^{(k)}$

    2. (b)

$Y\leftarrow Y+A^{(k)^{T}}\hat{U}-\hat{U}A^{(k)^{T}}$

  3. 3.

Return $\text{vec}(Y)$

The computational complexity can further be reduced if each of the matrices
$A^{(k)}$ is sparse with average density $\rho\in[0,1]$ (defined as the number
of nonzero entries divided by $N^{2}$) so computing the matrix-matrix products
requires in average only $4\rho N^{3}$.
The dominant computational complexity of each step of the Lanczos algorithm is
the matrix vector product Saad (2003) which we have shown can be computed in
$4\rho MN^{3}$ operations, far more cheaply than the $N^{4}$ operations if $S$
did not have its special structure. Due to the iterative nature of the Lanczos
algorithm, it is unknown a priori the number of iterations required to compute
the eigenvalue/eigenvector pairs. Nonetheless, unless the number of iterations
required is on the order of $N^{2}$ or larger, the Lanczos algorithm is more
efficient than constructing $S$ explicitly and finding the
eigenvalue/eigenvectors pairs using a standard tridiagonalization approach for
dense symmetric eigenvalue problems.
Let $\bm{v}_{k}\in\mathbb{R}^{NN}$, $k=1,2,\ldots,n_{ev}$, be the $n_{ev}$
eigenvectors found corresponding to the eigenvalues equal to zero of $S$, each
of which lies in the intersection of the nullspaces of $P^{(k)}$,
$k=1,2,\ldots,M$. To select a random vector in the intersection of the
nullspaces, create the vector
$\text{vec}(U)=\sum_{k=1}^{n_{ev}}c_{k}\bm{v}_{k}$ where we uniformly at
random select $c_{k}\in[-1,1]$, $k=1,2,\ldots,n_{ev}$, and scale them such
that $\sum_{k}c_{k}^{2}=1$. The resulting matrix $U$ satisfies all of the
commutation relations, $[A^{(k)},U]=O_{N}$, as does $U^{T}$ because each
$A^{(k)}$ is symmetric. With this fact, the symmetric matrix
$\frac{1}{2}(U+U^{T})$ also commutes,
$\left[A^{(k)},\frac{1}{2}(U+U^{T})\right]$, $k=1,2,\ldots,M$. Finally, to
find the matrix $P$ that simultaneously block diagonalizes each of the
$A^{(k)}$, $k=1,2,\ldots,M$, compute the eigenvectors of
$\frac{1}{2}(U+U^{T})$, and store them as the columns of $P$. The proof of the
correctness is extensive and beyond the scope of this paper but can be found
in Maehara and Murota (2011).
A related problem Maehara and Murota (2011) is to find a transformation $P$
that does not exactly simultaneously block diagonalize all of the matrices
$A^{(k)}$, $k=1,2,\ldots,M$, but rather results in matrices with off-diagonal
blocks with entries with magnitude of the order $\epsilon>0$. The process is
validated by Lemma 4.1 in Maehara and Murota (2011) and proceeds exactly as
before except now rather than finding the eigenvectors associated with
eigenvalues equal to zero, instead, the Lanczos method is used to find
eigenvalues of $S$ less than $\epsilon$ along with their eigenvectors. After
this, with the eigenvectors $\bm{v}_{k}$, $k=1,2,\ldots,n_{ev}$, the same
steps are taken to compute $U$, extract its symmetric part
$\frac{1}{2}(U+U^{T})$, and find the eigenvectors of the result.

We make available our code to compute the FSBD of a set of symmetric matrices
– see the Code Availability Section.

##  IV Application of the SBD technique to complete synchronization of
networks with different types of connections

The time evolution of a network of dynamical systems coupled through different
types of connections is described by the following set of equations:

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

where $\boldsymbol{x}_{i}(t)$ and $\boldsymbol{F}(\boldsymbol{x}_{i}(t))$
represent the $m$-dimensional state vector and dynamical function of the
system located at node $i$, respectively. The network nodes are coupled
through different coupling functions $\boldsymbol{H}^{(k)}$, $k=1,\cdots M$.
The network connectivity associated to each coupling function is described by
the adjacency matrix $A^{(k)}$, where $A_{ij}^{(k)}=A_{ji}^{(k)}>0$ if there
is a connection between nodes $i$ and $j$ and $A_{ij}^{(k)}=A_{ji}^{(k)}=0$
otherwise. The above set of equations can be rewritten,

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

where the Laplacian matrices $L^{(k)}$ have entries
$L_{ij}^{(k)}=A_{ij}^{(k)}>0$ for $j\neq i$ and $L_{ii}^{(k)}=-\sum_{j\neq
i}L_{ij}^{(k)}$, $k=1,...,M$. Note that all of the rows of the Laplacian
matrices $L^{(k)}$ sum to zero, $k=1,...,M$. The synchronization manifold
$\boldsymbol{x}_{1}=\boldsymbol{x}_{2}=\cdots=\boldsymbol{x}_{N}$ is an
invariant subspace for the set of Eqs. (7). The dynamics on this manifold,
which corresponds to complete synchronization,
$\boldsymbol{x}_{1}(t)=\boldsymbol{x}_{2}(t)=\cdots=\boldsymbol{x}_{N}(t)=\boldsymbol{x}_{s}(t)$
obeys the equation of an uncoupled system,

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

To investigate the stability of the complete synchronous state, we study the
dynamics of a small perturbation from the synchronous solution
$(\boldsymbol{x}_{i}(t)=\boldsymbol{x}_{s}(t)+\delta\boldsymbol{x}_{i}(t))$.
The synchronous state is stable if the perturbations approach $0$ for large
$t$. The linearized system of equations can be written,

$\delta\dot{\boldsymbol{x}}_{i}(t)=D\boldsymbol{F}(\boldsymbol{x}_{s}(t))\delta\boldsymbol{x}_{i}(t)+\sum_{k=1}^{M}\sum_{j=1}^{N}L_{ij}^{(k)}D\boldsymbol{H}^{(k)}(\boldsymbol{x}_{s}(t))\delta\boldsymbol{x}_{j}(t)\quad
i=1,\cdots,N$ (9)

By stacking together all the perturbations in an $mN$-dimensional vector
$\boldsymbol{z}=[\delta\boldsymbol{x}_{1}^{T},\delta\boldsymbol{x}_{2}^{T},\cdots,\delta\boldsymbol{x}_{N}^{T}]^{T}$,
the set of Eqs. (9) can be rewritten in vectorial form,

$\dot{\boldsymbol{z}}(t)=[I_{N}\otimes
D\boldsymbol{F}(\boldsymbol{x}_{s}(t))+\sum_{k=1}^{M}L^{(k)}\otimes
D\boldsymbol{H}^{(k)}(\boldsymbol{x}_{s}(t))]\boldsymbol{z}(t)$ (10)

where $\otimes$ is the Kronecker product.

One observation is that by construction the set of Laplacian matrices
$L^{(1)},L^{(2)},...,L^{(M)}$ all share one common eigenvector
$[1,1,...,1]/\sqrt{N}$, with associated eigenvalue $0$. It follows that we can
define an orthogonal transformation $\tilde{P}$ leading to a trivial
simultaneous block diagonalization (TSBD)

$\tilde{P}^{T}L^{(k)}\tilde{P}=0\oplus B_{\mbox{orth}},\quad k=1,...,M,$ (11)

where the block $B_{\mbox{orth}}$ is $(N-1)$-dimensional. Hence, there will be
a large block produced by the TSBD with dimension $N-1$. For example,
$\tilde{P}$ could be taken as the matrix whose columns are the eigenvectors of
any of the matrices $L^{(1)}$, $L^{(2)},...,L^{(M)}$ Sorrentino (2012). One
would hope that calculation of the FSBD for the set of matrices
$L^{(1)},L^{(2)},...,L^{(M)}$ leads to a finer block-diagonalization than the
TSBD.

Here, for simplicity and without loss of generality, we focus on the case of
$M=2$ different connection types, for which Eq. (10) becomes,

$\dot{\boldsymbol{z}}(t)=\left[I_{N}\otimes
D\boldsymbol{F}(\boldsymbol{x}_{s}(t))+L^{(1)}\otimes
D\boldsymbol{H}^{(1)}(\boldsymbol{x}_{s}(t))+L^{(2)}\otimes
D\boldsymbol{H}^{(2)}(\boldsymbol{x}_{s}(t))\right]\boldsymbol{z}(t)$ (12)

We attempt to break the stability problem (12) into a set of independent
lower-dimensional equations. To this end, we seek for a transformation that
leads to decoupling the set of Eqs. (12), by simultaneously block
diagonalizing $L^{(1)}$ and $L^{(2)}$. Special instances of this problem have
been studied in Ref. Sorrentino (2012) which obtained three different
conditions under which the problem with $mN$-dimension can be broken into a
set of $(N-1)$ problems of dimension $m$ each. Moreover, Ref. Irving and
Sorrentino (2012) has introduced the general framework in which the SBD
technique is applied to network synchronization.

We now compute $P=\mathcal{SBD}(L^{(1)},L^{(2)})$ and rewrite Eq. (12) as
follows:

$\dot{\boldsymbol{\eta}}(t)=\left[I_{N}\otimes
D\boldsymbol{F}(\boldsymbol{x}_{s}(t))+(P^{T}L^{(1)}P)\otimes
D\boldsymbol{H}^{(1)}(\boldsymbol{x}_{s}(t))+(P^{T}L^{(2)}P)\otimes
D\boldsymbol{H}^{(2)}(\boldsymbol{x}_{s}(t))\right]\boldsymbol{\eta}(t),$ (13)

where $\boldsymbol{\eta}(t)=(P\otimes I_{m})\boldsymbol{z}(t)$.

As stated before, $P^{T}L^{(i)}P=\oplus_{j=1}^{n}B_{j}^{(i)}$, where all the
matrices $B_{j}$ have the same block-diagonal form. Therefore, Eqs. (13) can
be decoupled as follows,

$\dot{\boldsymbol{\eta}}_{i}(t)=\left[I_{D_{i}}\otimes
D\boldsymbol{F}(\boldsymbol{x}_{s}(t))+B_{i}^{(1)}\otimes
D\boldsymbol{H}^{(1)}(\boldsymbol{x}_{s}(t))+B_{i}^{(2)}\otimes
D\boldsymbol{H}^{(2)}(\boldsymbol{x}_{s}(t))\right]\boldsymbol{\eta}_{i}(t),$
(14)

where $D_{i}$ is the block-dimension of $B_{i}$, $\sum_{i}D_{i}=N$.

Recall that for a given value of $i$ ($i=1$) we obtain scalar blocks
$B^{(1)}_{1}=B^{(2)}_{1}=0$, which are associated with a perturbation parallel
to the synchronization manifold. Therefore, to analyze the stability of the
synchronous solution, we only need to deal with the other $i>1$ transverse
blocks (14).

###  IV.1 Performance of the SBD technique applied to complete synchronization

To examine the extent of the reduction provided by the SBD method, the
following index is introduced

$0\leq I_{d}=1-\frac{Z-1}{N-2}\leq 1$ (15)

where $Z$ is the maximum block dimension in both $P^{T}L^{1}P$ and
$P^{T}L^{2}P$. The best possible performance of the SBD is achieved for
$I_{d}=1$, corresponding to all the blocks having dimension $Z=1$. On the
other hand, if the maximum block dimension $Z=N-1$, the index $I_{d}=0$, which
corresponds to the same reduction achievable with the TSBD.

In the following, we examine the performance of the SBD method in reducing the
problem of complete synchronization for three different networks: (i) Erdős-
Rényi (ER) random networksErdős and Rényi (1960), (ii) Watts-Strogatz small-
world (WS) networksWatts and Strogatz (1998), and (iii) scale-free networks
Barabási and Bonabeau (2003) generated by the configuration model Molloy and
Reed (1995). For ER networks, the adjacency matrices $A^{(1)}$ and $A^{(2)}$
are taken to describe the connectivity of two ER networks with connection
probabilities $p_{1}$ and $p_{2}$, respectively. After computing the SBD for
the two Laplacian matrices $L^{(1)}$ and $L^{(2)}$, the index $I_{d}$ is
calculated to investigate the performance of the SBD algorithm. Fig. 1 a. b.
and c. show the index $I_{d}$ versus the probabilities $p_{1}$ and $p_{2}$ for
the network with $N=10,\>20$ and $50$ nodes, respectively. The connectivity of
two WS small-world networks are shown with two adjacency matrices $A^{(1)}$
and $A^{(2)}$ in which the connection probabilities $p_{1}$ and $p_{2}$ are
chosen in the interval $[0.01\quad 0.1]$ to fulfill the WS small-world
properties. Fig. 1 e. and f. show the index $I_{d}$ versus the probabilities
$p_{1}$ and $p_{2}$ for the WS small-world network with $N=20$ and $N=50$
nodes, respectively. Fig. 1 d shows the index $I_{d}$ versus the exponent of
power-law distribution $\gamma_{1}$ and $\gamma_{2}$ of scale-free networks
with $N=50$ nodes where the minimum degree of each node is set to $k=3$ in
order to have a connected network.

Figure 1: a-c) The index $I_{d}$ for two Erdős-Rényi topologies is plotted vs
the connection probabilities $p_{1}$ and $p_{2}$. In a) $N=10$ nodes, in b)
$N=20$ nodes, and in c) $N=50$ nodes. d) The index $I_{d}$ for two scale-free
networks with $N=50$ nodes is plotted vs the power-law distributions
$\gamma_{1}$ and $\gamma_{2}$. e and f) The index $I_{d}$ for two Watts-
Strogatz small-world topologies is plotted vs the connection probabilities
$p_{1}$ and $p_{2}$ in networks with $N=20$ and $N=50$ nodes. Different values
of the index $I_{d}$ are shown as variations in the color spectrum from dark
blue ($I_{d}=0$) to dark red ($I_{d}=1$).

The size of the networks in Panel a) is $N=10$ nodes, in Panel b and e) is
$N=20$ nodes, and in Panel c, d and f) is $N=50$ nodes. Different values of
the index $I_{d}$ are shown as variations in the color spectrum from dark blue
($I_{d}=0$) to dark red ($I_{d}=1$). As can be seen, most of the parameter
space is colored dark blue. In panel a-c) and e-f) increasing the
dimensionality of the network from $N=10$ to $N=50$ corresponds to increasing
the extension of the dark blue area. Also, Panel d) shows that different
values of the exponent of power-law distribution of two scale-free networks
with $50$ nodes between the interval of $[2\quad 5]$ results in a large dark
blue area. All of our results show that the SBD technique fails when applied
to two random networks since for most choices of $p_{1}$ and $p_{2}$ or
$\gamma_{1}$ and $\gamma_{2}$, we obtain $I_{d}=0$.
To be more precise, there are some areas in Fig. 1 Panel a-c) and e-f)
corresponding to $I_{d}>0$. These areas are associated with very small and
very large connection probabilities $p_{i}$, i.e., with the limit of either
very sparse ER or small-world network or very dense ER network. For instance,
consider the ER networks with two limits in which $p_{i}=0$ and $p_{i}=1$. For
$p_{i}=0$, the corresponding adjacency matrix is a zero matrix. Then, the
stability problem is reduced to diagonalization of just the other adjacency
matrix. For the case $p_{i}=1$, the corresponding ER network is a fully-
connected network which corresponds to one of the dimensionality reduction
studied in Ref. Sorrentino (2012).

##  V Application of the SBD technique to cluster synchronization

The problem of investigating the stability of cluster synchronous solutions in
networks has attracted much attention in the last few years. A general
equation for a network of coupled dynamical systems is the following,

$\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$ (16)

where the network connectivity is described by the adjacency matrix $A$, where
$A_{ij}=A_{ji}>0$ if there is a connection between nodes $i$ and $j$ and
$A_{ij}=A_{ji}=0$ otherwise. The function $\boldsymbol{H}$ is the node-to-node
coupling function.

The nodes of the network can be partitioned in a set of $C$ equitable clusters
or balanced colors $\mathcal{C}_{1},\mathcal{C}_{2},...,\mathcal{C}_{C}$,
where $N_{i}$ is the number of nodes in cluster $\mathcal{C}_{i}$ and
$\sum_{i=1}^{C}N_{i}=N$ Sachs (1966); Schaub  _et al._ (2016a). All the nodes
in the same equitable cluster receive the same number of connections from each
one of the clusters Egerstedt  _et al._ (2012). Among several possible
equitable partitions of the network, there is one corresponding to the minimum
number of clusters, which we will refer to as the minimum balanced coloring.
For any adjacency matrix $A$, the algorithm described by Belykh and Hasler
Belykh and Hasler (2011) outputs the minimum balanced coloring very
efficiently. Information about the minimum balanced coloring is contained in
the $N\times C$ indicator matrix $O=\\{O_{ij}\\}$ where $O_{ij}$ is equal to 1
if node $i$ is in cluster $C_{j}$ and is 0 otherwise.

Similarly to the case of complete synchronization considered earlier, given an
equitable partition of the network nodes, we can define an invariant subspace
for the set of Eqs. (16), 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{x}^{s}_{1}(t),\boldsymbol{x}^{s}_{2}(t),\cdots,\boldsymbol{x}^{s}_{C}(t)\\}$,
where $\boldsymbol{x}^{s}_{1}(t)$ is the synchronous solution for nodes in
cluster $\mathcal{C}_{1}$, $\boldsymbol{x}^{s}_{2}(t)$ is the synchronous
solution for nodes in cluster $\mathcal{C}_{2}$, and so on.
We can then define the $C\times C$ quotient matrix $Q$ such that for each pair
of clusters $\mathcal{C}_{u}$ and $\mathcal{C}_{v}$,

$Q_{uv}=\sum_{j\in\mathcal{C}_{v}}A_{ij}\quad i\in\mathcal{C}_{u}\quad
u,v=1,2,\cdots,C$ (17)

All the nodes belonging to the same cluster can synchronize on the quotient
network time evolution $(\boldsymbol{x}_{u}^{s}(t))$, $u=1,...,C$

$\dot{\boldsymbol{x}_{u}^{s}}(t)=\boldsymbol{F}(\boldsymbol{x}_{u}^{s}(t))+\sum_{v=1}^{C}Q_{uv}\boldsymbol{H}(\boldsymbol{x}_{v}^{s}(t)),$
(18)

$u=1,...,C$. The question we are interested in is whether the cluster
synchronous solution corresponding to the minimum balanced coloring is stable
or unstable.

Stability of the cluster synchronous solution depends on the $mN$-dimensional
equation,

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

where $E_{c}$ is a diagonal matrix such that $(E_{c})_{ii}=1$ if node $i$
belongs to cluster $c$ and $(E_{c})_{ii}=0$ otherwrise.

We note that by left-multiplying Eq. (19) by the matrix $\tilde{O}\otimes
I_{m}$ where $\tilde{O}=(O^{T}O)^{-1}O^{T}$ we obtain the dynamics of the
perturbation parallel to the synchronization manifold Schaub  _et al._
(2016b); Siddique  _et al._ (2018); Klickstein, Pecora, and Sorrentino (2019).

Similarly to Sec. IV, we would like to reduce the stablility problem into a
set of independent lower-dimensional equations instead of dealing with the
high-dimensional problem, Eq. (19). Ref. Pecora  _et al._ (2014) has proposed
a dimensionality reduction approach based on group theory for the case of
orbital clusters and shown that the IRR of the symmetry group can be used to
block-diagonalize the set of Eq. (19). Ref. Zhang and Motter (2020) has
applied the SBD method to characterize stability of any cluster
synchronization pattern. An important question is whether the symmetry-
independent approach of Zhang and Motter (2020) may lead to a dimensionality
reduction of the stability analysis in the broader class of networks
Klickstein and Sorrentino (2018a) that have equitable clusters that are not
merely the result of symmetries Siddique  _et al._ (2018). Next we show that
this is not the case.

Following Zhang and Motter (2020) we compute
$P=\mathcal{SBD}(A,E_{1},E_{2},...,E_{C})$. By applying $P$ to Eq. (19), we
obtain,

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

where $\boldsymbol{\eta}(t)=P^{T}\otimes I_{m}\boldsymbol{z}(t)$. Note that
because both matrices $P^{T}AP$ and $P^{T}E_{c}P$ have the same block diagonal
structure, so does $P^{T}AE_{c}P$, which becomes apparent by rewriting
$P^{T}AE_{c}P=(P^{T}AP)(P^{T}E_{c}P)$. Therefore, (20) can be decomposed into
smaller blocks,

$\dot{\boldsymbol{\eta}}_{i}(t)=\left[\sum_{c=1}^{C}(J_{i})_{c}\otimes
D\boldsymbol{F}(\boldsymbol{x}^{s}_{c}(t))+\sum_{c=1}^{C}(B_{i})_{c}\otimes
D\boldsymbol{H}^{2}(\boldsymbol{x}^{s}_{c}(t))\right]{\boldsymbol{\eta}}_{i}(t)$
(21)

where $(J_{i})_{c}$ and $(B_{i})_{c}$ are blocks of the same dimensions
derived from the transformations
$P^{T}\sum_{c=1}^{C}E_{c}P=\otimes_{j=1}^{n}(J_{j})_{c}$ and
$P^{T}\sum_{c=1}^{C}AE_{c}P=\otimes_{j=1}^{n}(B_{j})_{c}$, respectively.

###  V.1 Generating networks with assigned equitable partition

In order to study the performance of the SBD reduction in the case of cluster
synchronization, we need a method to generate a random symmetric network with
an assigned equitable partition. This can be done by using the algorithm
described below.

First, assign the number of nodes in each of $C$ clusters,
$N_{1},N_{2},...,N_{C}$, where in order to enforce a trivial pattern of
connectivity we pick $N_{1},N_{2},...,N_{C}$, so that no two such numbers are
coprime, i.e. $\mbox{mcd}(N_{i},N_{j})>1$, $i=1,...,C$, $j\neq i$. Second, we
need to determine the relative indegree $d_{ij}$ of nodes in cluster $i$ from
nodes in cluster $j$. Due to the assumption that the network is symmetric, the
following condition needs to be satisfied

$N_{i}d_{ij}=N_{j}d_{ji}.$ (22)

A trivial solution is $d_{ij}=N_{j}$ and $d_{ji}=N_{i}$, which corresponds to
a trivial pattern of connectivity in which each node in cluster $i$ is coupled
to all the nodes in cluster $j$ and vice versa. By the assumption that $N_{i}$
and $N_{j}$ are not coprimes, it follows that we can always choose nontrivial
values of $d_{ij}$ and $d_{ji}$, namely,

$d_{ij}=\frac{N_{j}}{\alpha}\qquad d_{ji}=\frac{N_{i}}{\alpha},$ (23)

where $\alpha=\mbox{mcd}(N_{i},N_{j})>1$. Then, for each pair of clusters, we
can randomly connect the nodes in cluster $i$ and cluster $j$ with
$\frac{N_{i}N_{j}}{\alpha}$ bidirectional links. The intra-connectivity of
each cluster is determined by first assigning the intra-degree $d_{i}$ of all
nodes in cluster $i$ for $i=1,2,...,N_{C}$. This should be chosen such that
$N_{i}d_{i}$ is an even number and $d_{i}<N_{i}$. The latter condition will
assure that nodes $n_{i}$ and $n_{j}$ for $i\neq j$ share no more than one
connection.

The algorithm provided here generates a network with an assigned equitable
partition as opposed to the algorithms to generate networks with assigned
orbital partition presented in Klickstein and Sorrentino (2018a, b).

###  V.2 Performance of the SBD technique applied to cluster synchronization

For the case of cluster synchronization, we can also define a transformation
corresponding to the trivial simultaneous block diagonalization (TSBD.) This
corresponds to the transformation that separates the perturbation parallel to
the synchronization manifold from the perturbation transverse to the
synchronization manifold. By choosing $\tilde{P}=\oplus_{i=1}^{C}G_{i}$ where
$G_{i}$ is an orthogonal matrix of dimension $N_{i}$ with one of its columns
equal to $\textbf{1}_{N_{i}}/\sqrt{N_{i}}$ we obtain the trivial simultaneous
block diagonalization

$P^{T}AP=B_{\mbox{par}}\oplus B_{\mbox{orth}},$ (24)

where the block $B_{\mbox{par}}$ is $C$-dimensional and the block
$B_{\mbox{orth}}$ is $(N-C)$-dimensional. Hence, the largest block produced by
the TSBD will have dimension $L=\max{(C,N-C)}$. Thus for the case of cluster
synchronization, we define the performance index,

$I_{cs}=\frac{L-Z}{L}$ (25)

Where $Z$ is the largest block dimension resulting from calculation of the
FSBD for the set of matrices $\\{A,E_{1},E_{2},...,E_{C}\\}$ where $A$ is the
adjacency matrix and $E_{1},E_{2},...,E_{C}$ are the matrices containing
cluster information. Again the index compares the performance of the FSBD with
that of the TSBD. An index $I_{cs}=0$ indicates that the reduction achieved by
the FSBD is the same as that of the TSBD and is therefore regarded as a
failure.

Next we consider a numerical example for a random symmetric network with $C=4$
clusters generated using the algorithm previously described. Figure 2 shows
the network that is to be examined to measure the performance index for the
SBD algorithm, with nodes color coded according to the equitable cluster they
belong to.

Figure 2: Randomly constructed symmetric network with $C=4$ equitable clusters
and 34 nodes. The clusters are identified as follows: $\mathcal{C}_{1}$ (green
nodes) with $N_{1}=4$, $\mathcal{C}_{2}$ (red nodes) with $N_{2}=8$,
$\mathcal{C}_{3}$ (black nodes) with $N_{3}=12$, and $\mathcal{C}_{4}$
(magenta nodes) with $N_{4}=10$. The arbitrarily chosen intra-degrees for this
network are $d_{1}=2$, $d_{2}=3$, $d_{3}=4$, and $d_{4}=6$.

Figure 3 shows block-diagonalized matrices obtained by application of the FSBD
procedure to the network in Fig. 2.

Figure 3: Block-diagonalized matrices for the randomly constructed symmetric
network shown in Fig 2. By applying FSBD to the set
$\\{A,E_{1},E_{2},E_{3},E_{4}\\}$ the new set
$\\{P^{T}AP,P^{T}E_{1}P,P^{T}E_{2}P,P^{T}E_{3}P,P^{T}E_{4}P\\}$ is obtained.
Each blue dot represents a non-zero entry in the new set.

In order to better visualize the block decomposition, we construct the matrix
$\Omega$ as the sum of absolute values of the matrices
$\\{P^{T}AP,P^{T}E_{1}P,P^{T}E_{2}P,P^{T}E_{3}P,P^{T}E_{4}P\\}$

$\Omega=|P^{T}AP|+|P^{T}E_{1}P|+|P^{T}E_{2}P|+|P^{T}E_{3}P|+|P^{T}E_{4}P|$
(26)

where the symbol $|\cdot|$ here indicates the entry-wise absolute value of a
matrix. A representation of the matrix $\Omega$ is shown in Fig. 4, which
evidences two blocks: one $4$-dimensional block, and one $30$-dimensional
block. For this example, the calculated performance index is $0$ with $L=30$
and $Z=30$. We have obtained similar results for all the other instances of
random graphs with assigned equitable partition (algorithm of Sec. VA) we have
tested.

Figure 4: Representation of the matrix $\Omega$ for the network shown in Fig.
2. The red dots represent non-zero entries. Using equation (25), the
calculated performance index is $I_{cs}=0$ with $L=30$ and $Z=30$ for this
example.

##  VI Conclusions

The techniques for simultaneous block diagonalization of matrices have been
developed by Maehara, Murota _et al_ in a number of seminal papers Maehara and
Murota (2010a, b); Murota  _et al._ (2010); Maehara and Murota (2011). These
techniques were originally applied to problems in the areas of semidefinite
programming and signal processing (independent component analysis), see e.g.
Maehara and Murota (2010a). The first application of these techniques to
network synchronization was presented in a 2012 paper Irving and Sorrentino
(2012). Only recently they have been applied to the problem of cluster
synchronization of networks Zhang and Motter (2020).

We are highly indebted to the mathematicians who have developed the algebraic
theory of simultaneous block diagonalization of matrices. This is relevant to
many problems in the applied sciences where one is looking for modal
decompositions but such decompositions may not be obvious. The application of
these techniques to the problem of network synchronization is important as it
allows to define the extent to which the synchronization stability problem can
be reduced in realistic situations that deviate from the original assumptions
of nodes all of the same type and connections all of the same type Pecora and
Carroll (1998a). We have seen in this paper that unfortunately in generic
situations (random networks) the obtained reduction is modest and comparable
to that achievable with a trivial transformation. Even if that is the case, it
is important to know the extent of the attainable reduction and that no
further decomposition of the problem is possible. With this paper we believe
we have set the expectations straight about the reduction that is
realistically achievable from application of SBD to the study of complete and
cluster synchronization of generic (random) graphs. Overall, this does not
diminish our enthusiasm for these techniques, which can provide exceptional
insight into many problems of interest in physics and engineering, including
network synchronization. Moreover, both Refs. Irving and Sorrentino (2012) and
Zhang and Motter (2020) have shown that the reduction produced by the SBD
technique can be substantial for specific network realizations, which provides
an important and valuable result when one has the ability to appropriately
design the networks connectivity.

##  VII Code Availability

Code to compute the simultaneous block diagonalizations for the examples shown
in this paper can be accessed from its  Github repo.

## References

  * Bai  _et al._ (2009) Bai, Y., de Klerk, E., Pasechnik, D., and Sotirov, R., “Exploiting group symmetry in truss topology optimization,” Optimization and Engineering 10, 331–349 (2009).
  * Barabási and Bonabeau (2003) Barabási, A.-L. and Bonabeau, E., “Scale-free networks,” Scientific american 288, 60–69 (2003).
  * 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).
  * Bernstein (2009) Bernstein, D. S.,  _Matrix mathematics_ (Princeton university press, 2009).
  * De Klerk and Sotirov (2010) De Klerk, E. and Sotirov, R., “Exploiting group symmetry in semidefinite programming relaxations of the quadratic assignment problem,” Mathematical Programming 122, 225–246 (2010).
  * Egerstedt  _et al._ (2012) Egerstedt, M., Martini, S., Cao, M., Camlibel, K., and Bicchi, A., “Interacting with networks: How does structure relate to controllability in single-leader, consensus networks?” IEEE control systems magazine 32, 66–73 (2012).
  * Erdős and Rényi (1960) Erdős, P. and Rényi, A., “On the evolution of random graphs,” Publ. Math. Inst. Hung. Acad. Sci 5, 17–60 (1960).
  * Gatermann and Parrilo (2004) Gatermann, K. and Parrilo, P. A., “Symmetry groups, semidefinite programs, and sums of squares,” Journal of Pure and Applied Algebra 192, 95–128 (2004).
  * 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 a hypernetwork of coupled dynamical systems,” Phys. Rev. E 86, 056102 (2012).
  * Kanno  _et al._ (2001) Kanno, Y., Ohsaki, M., Murota, K., and Katoh, N., “Group symmetry in interior-point methods for semidefinite program,” Optimization and Engineering 2, 293–320 (2001).
  * 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).
  * Klickstein and Sorrentino (2018a) Klickstein, I. and Sorrentino, F., “Generating graphs with symmetry,” IEEE Transactions on Network Science and Engineering 6, 836–843 (2018a).
  * Klickstein and Sorrentino (2018b) Klickstein, I. and Sorrentino, F., “Generating symmetric graphs,” Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 121102 (2018b).
  * Kojima, Kojima, and Hara (1997) Kojima, M., Kojima, S., and Hara, S., “Linear algebra for semidefinite programming,” 1004, 1–23 (1997).
  * 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).
  * Miller (1973) Miller, W.,  _Symmetry groups and their applications_ (Academic Press, 1973).
  * Molloy and Reed (1995) Molloy, M. and Reed, B., “A critical point for random graphs with a given degree sequence,” Random Structures and Algorithms , 161–179 (1995).
  * Murota  _et al._ (2010) Murota, K., Kanno, Y., Kojima, M., and Kojima, S., “A numerical algorithm for block-diagonal decomposition of matrix $\ast$-algebras with application to semidefinite programming,” Japan Journal of Industrial and Applied Mathematics 27, 125–160 (2010).
  * Pecora and Carroll (1998a) Pecora, L. and Carroll, T., “Master stability functions for synchronized coupled systems,” Phys. Rev. Lett. 80, 2109–2112 (1998a).
  * Pecora and Carroll (1998b) Pecora, L. M. and Carroll, T. L., “Master stability functions for synchronized coupled systems,” Physical review letters 80, 2109 (1998b).
  * 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).
  * Riener  _et al._ (2013) Riener, C., Theobald, T., Andrén, L. J., and Lasserre, J. B., “Exploiting symmetries in sdp-relaxations for polynomial optimization,” Mathematics of Operations Research 38, 122–141 (2013).
  * Saad (2003) Saad, Y.,  _Iterative methods for sparse linear systems_ (SIAM, 2003).
  * Sachs (1966) Sachs, H., “Über teiler, faktoren und charakteristische polynome von graphen,” Teil I. Wiss. Z. TH Ilmenau 12, 7–12 (1966).
  * Schaub  _et al._ (2016a) 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 (2016a).
  * Schaub  _et al._ (2016b) 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 (2016b).
  * Serre (1977) Serre, J.-P.,  _Linear representations of finite groups_ , Vol. 42 (Springer, 1977).
  * 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 (2012) Sorrentino, F., “Synchronization of hypernetworks of coupled dynamical systems,” New J. Phys. 14, 033035 (2012).
  * Uhlig (1973) Uhlig, F., “Simultaneous block diagonalization of two real symmetric matrices,” Linear Algebra and its Applications 7, 281–289 (1973).
  * Watts and Strogatz (1998) Watts, D. J. and Strogatz, S. H., “Collective dynamics of ‘small-world’networks,” nature 393, 440–442 (1998).
  * Wedderburn (1934) Wedderburn, J. H. M.,  _Lectures on matrices_ , Vol. 17 (American Mathematical Soc., 1934).
  * Zhang and Motter (2020) Zhang, Y. and Motter, A. E., “Symmetry-independent stability analysis of synchronization patterns,” SIAM Review 62, 817–836 (2020).