A polynomial eigenvalue approach for multiplex networks

Guilherme Ferraz de Arruda; Emanuele Cozzo; Francisco A Rodrigues; Yamir Moreno

doi:10.1088/1367-2630/aadf9f

1. Introduction

Complex network theory has become one of the main tools for the analysis of complex systems, allowing the representation of a wide range of systems composed by interacting discrete elements [1]. However, real-life systems are also organized in layers, which represent different channels of interaction. In order to incorporate these characteristics, one should work with multilayer networks, which allows for a proper representation of multiplex and interconnected systems [2–4]. The introduction of this extra level of complexity also imposes new challenges on the analysis of its structural and dynamical properties. Furthermore, a key element on the analysis of networks is their spectral properties [5]. In fact, they play an important role in explaining the connection between structure and dynamics. For instance, in epidemic spreading the critical point below which the infection prevalence is null is predicted to be the inverse of the leading eigenvalue of the adjacency matrix, in both, single [6] and multiplex networks [7]. Additionally, its nature also seems to be connected to those properties [7–9]. Although the literature about the spectra of single-layer networks is well developed [5], the theory of spectral properties of multiplex networks is still in its infancy. This motivates us to propose a different formalism aimed at filling this gap.

In this paper, we will consider the matrix representation of multiplex networks, constraining ourselves to finite matrices. First of all, we are interested in weighting differently inter and intra-layer edges. This implies that those matrices will also be a function of the inter-layer coupling parameter, here called p. Consequently, the associated eigensystem will be a function of that same parameter. Additionally, the matrix approach is especially interesting in this context since it allows us to directly use linear algebra and spectral graph results already available.

When varying the coupling parameter, a multiplex system might present different structural phases, which are characterized in terms of eigenvalue crossings and eigengaps and are intuitively defined as: (i) decoupled phase, for small values of p, where the layers are virtually decoupled and act by themselves, with a negligible interaction between layers, (ii) multiplex/multilayer phase, where the system is coupled and the intra-layer edges play an important role and (iii) the aggregate network phase, where the system behaves as the superposition of all layers. It is clear that a good understanding of the eigenvalues' behavior might be useful since we could move our system into different structural regimes, aiming at different goals such as improved robustness, better performance regarding diffusion or spreading, among many other possible applications.

The different structural phases above are related through the interlacing properties of quotient graphs [10, 11]. More specifically, in [10] the authors showed that the spectra of different scales of a multiplex network (aggregated network, the network of layers and individual layers) characterize the three phases. In practical terms, the interlacing provides us bounds for the spectra [10, 11], but also emphasizes that the different scales are intrinsically connected. Indeed, it is impossible to tune the leading eigenvalue of the network of layers without also increasing the leading eigenvalue of the whole multiplex. Furthermore, in [12] the authors characterized multiple topological scales using the supra-Laplacian matrix. More specifically, they analyzed eigengaps to characterize them.

Following a different approach, in [13], the author evaluated the normalized Laplacian matrix, which, in fact, shares the same set of eigenvalues of the probability transition matrix (see section 4.3), and proposed a similar classification. However, it is worth mentioning that in [13], a different nomenclature was used and a fourth phase was defined. Namely, the proposed structural regimes were: (i) bipartite phase, (ii) decoupled phase, (iii) indistinguishable, where the author argues that the system is topologically and dynamically indistinguishable[13], and (iv) a mixed phase, (called BD in [13]—bipartite and decoupled phases) where the layers are structurally and dynamically distinguishable. Although here we do not make a distinction between the regimes (iii) and (iv) and consider both as multiplex regimes, we acknowledge the differences pointed out in [13]. It is also noteworthy that [13] considered structural correlations for the analysis, which is a key ingredient for the reported results. On the other hand, here we focus on uncorrelated networks and a multiplex structure.

The paper is organized as follows. In section 2, we present the polynomial eigenvalue formalism, giving its general definitions and properties in section 2.1. Next, we formalize the 2-Layer problem into a quadratic eigenvalue problem in section 3, analytically exploring its behavior as a function of the coupling parameter. In section 3.1, obtaining some bounds, in section 3.2, and discussing the simplified symmetric problem in section 3.3. We present our main applications in section 4, where we explore the supra-Laplacian matrix, in section 4.1, the supra-adjacency matrix in section 4.2 and the probability transition matrix in section 4.3. To round off this paper, we discuss the physical consequences of our findings, summarize our main results and perspectives in section 5.

2. Polynomial eigenvalue problem

In this section, we formally define the polynomial eigenvalue problem and present some of its fundamental properties. The aim of this section is to generically define our main mathematical object, establishing its basic properties. Thus, this will allow us to properly study the matrices associated with 2-Layer multiplex networks, which will naturally appear as a consequence of a simple manipulation of a linear system describing the network. From this simple approach, we expect to provide a different perspective on the spectral properties of multiplex networks.

2.1. General definition and properties

A matrix polynomial of order l is a matrix-valued function of a complex variable of the form [14]

$\begin{eqnarray}&&{\bf{R}}(\lambda )=\displaystyle \sum _{i=0}^{l}{{\bf{M}}}_{i}{\lambda }^{i},\end{eqnarray} \tag{ 1 }$

where ${{\bf{M}}}_{0},{{\bf{M}}}_{1},\,\ldots ,\,{{\bf{M}}}_{l}$ are n × n matrices and they are said coefficient matrices. If ${{\bf{M}}}_{l}={\bf{I}}$ the identity matrix the polynomial ${\bf{R}}$ is said to be monic. The eigenvalues of ${\bf{R}}$ are the solution to the characteristic equation

$\begin{eqnarray}&&\det ({\bf{R}}(\lambda ))=0.\end{eqnarray} \tag{ 2 }$

Therefore, right and left eigenvectors are defined as

$\begin{eqnarray}&&{\bf{R}}(\lambda )x=0,\end{eqnarray} \tag{ 3 }$

$\begin{eqnarray}&&{y}^{T}{\bf{R}}(\lambda )=0,\end{eqnarray} \tag{ 4 }$

where x and y are the right and left eigenvectors associated to the eigenvalue λ. It reduces to the standard eigenvalue problem when ${\bf{R}}(l)$ is a monic matrix polynomial with l = 1 and ${{\bf{M}}}_{0}=-{\bf{A}}$ , for any matrix ${\bf{A}}$ . A generalization of the Jordan form theory to a general matrix polynomial is possible and is briefly presented in the appendix.

A particular case of general interest is the quadratic eigenvalue problem (QEP), which is directly related to the 2-Layer case of our interest in this work. A quadratic matrix polynomial can be written as [14–16]

$\begin{eqnarray}&&{\bf{Q}}(\lambda )={\bf{A}}{\lambda }^{2}+{\bf{B}}\lambda +{\bf{C}}.\end{eqnarray} \tag{ 5 }$

Besides, without loss of generality, we assume that the eigenvectors are unitary. However, observe that if the right eigenvector x is unitary, thus, the left eigenvector y is not. Therefore, here we assume that x is unitary, which simplifies the equations.

Furthermore, if ${\bf{A}},{\bf{B}}$ and ${\bf{C}}$ are Hermitian a special class of problems, called hyperbolic quadratic eigenvalue problem (HQEP) [16], is obtained. Unfortunately, most of the problems in our context does not fall in this class.

3. A general 2-Layer case: a block matricial problem

The general form of any matrix associated to a multiplex network composed by two layers (supra-adjacency, supra-Laplacian and transition matrices for example) can be written as a block matrix. The resulting eigenvalue problem is the following

$\begin{eqnarray}\left[\begin{array}{cc}{{\bf{M}}}_{11} & {{\bf{M}}}_{12}\\ {{\bf{M}}}_{21} & {{\bf{M}}}_{22}\end{array}\right]\left[\begin{array}{c}{v}_{1}\\ {v}_{2}\end{array}\right]=\lambda \left[\begin{array}{c}{v}_{1}\\ {v}_{2}\end{array}\right],\end{eqnarray} \tag{ 6 }$

where ${{\bf{M}}}_{12}={{\bf{M}}}_{21}^{T}$ . Recognizing it as a system of equations and isolating v₁ on the second row, we have

$\begin{eqnarray}&&{v}_{1}=-{{\bf{M}}}_{21}^{-1}({{\bf{M}}}_{22}-\lambda {\bf{I}}){v}_{2}.\end{eqnarray} \tag{ 7 }$

Thus, inserting it in the first row we have

$\begin{eqnarray}&&[({{\bf{M}}}_{11}-\lambda {\bf{I}}){{\bf{M}}}_{21}^{-1}({{\bf{M}}}_{22}-\lambda {\bf{I}})-{{\bf{M}}}_{12}]{v}_{2}=0.\end{eqnarray} \tag{ 8 }$

This expression defines a QEP, whose coefficient matrices are

$\begin{eqnarray}&&{\bf{A}}={{\bf{M}}}_{21}^{-1},\end{eqnarray} \tag{ 9 }$

$\begin{eqnarray}&&{\bf{B}}=-({{\bf{M}}}_{11}{{\bf{M}}}_{21}^{-1}+{{\bf{M}}}_{21}^{-1}{{\bf{M}}}_{22}),\end{eqnarray} \tag{ 10 }$

$\begin{eqnarray}&&{\bf{C}}={{\bf{M}}}_{11}{{\bf{M}}}_{21}^{-1}{{\bf{M}}}_{22}-{{\bf{M}}}_{12},\end{eqnarray} \tag{ 11 }$

which poses a restriction on the inter-layer coupling matrix ${{\bf{M}}}_{12}$ , i.e., it must be invertible. Additionally, in our context, exchanging ${{\bf{M}}}_{11}$ and ${{\bf{M}}}_{22}$ we do not change the system, nor their solutions. Note that this operation is equivalent to relabeling the layers. However, denoting ${\bf{Q}}(\lambda )$ the polynomial of the first, thus, for the second case the polynomial is simply ${\bf{Q}}{(\lambda )}^{T}$ . In this way we also establish a relationship between the right and left eigenvectors and these two possible configurations. Formally, it implies that x = v₂ and y = v₁. As usual, we consider coupling matrices that are functions of a coupling parameter, p, i.e., ${{\bf{M}}}_{{ij}}={{\bf{M}}}_{{ij}}(p)$ , for $i\ne j$ . In fact, throughout this paper we explore how the spectral properties of our network evolve as we change such coupling parameter. As a constraint, we should mention that we only consider finite matrices.

Furthermore, note that ${\bf{B}}$ in equation (9) is intimately related to the aggregated and the loop-less aggregated networks of the original multiplex network (for more, see [10] or section 2.3.2 of [17]). More specifically, if the coupling matrices ${{\bf{M}}}_{12}$ and ${{\bf{M}}}_{21}$ are the identity matrix (or proportional to this matrix), thus ${\bf{B}}=-({{\bf{M}}}_{11}+{{\bf{M}}}_{22})$ , which is proportional to the loop-less aggregated network. Besides, note that a network of layers in the two-layer multiplex is a simple line graph with two nodes.

3.1. Spectral analysis as a function of p

So far, we made as less constraints as possible. Now we restrict ourselves to multiplex networks, i.e. diagonal coupling matrices, and, additionally, we assume a linear function of the parameter p > 0⁵ , ${{\bf{M}}}_{12}=p{\bf{D}}$ , where ${\bf{D}}$ is a diagonal invertible matrix (this constraint will be relaxed later). Then, we can describe each eigenvalue of ${\bf{Q}}(\lambda )$ by a scalar equation defined as the product of ${\bf{Q}}(\lambda )$ , by its left and right eigenvectors. Formally given as

$\begin{eqnarray}&&{y}^{T}{\bf{Q}}(\lambda )x=a({y}^{T},x){\lambda }^{2}+b({y}^{T},x)\lambda +c({y}^{T},x)=0,\end{eqnarray} \tag{ 12 }$

where $a({y}^{T},x)={y}^{T}{\bf{A}}x,b({y}^{T},x)={y}^{T}{\bf{B}}x$ and $c({y}^{T},x)={y}^{T}{\bf{C}}x$ . The solution of this equation is given by

$\begin{eqnarray}&&{\lambda }^{\pm }(x)=\displaystyle \frac{-b({y}^{T},x)\pm \sqrt{{\rm{\Delta }}({y}^{T},x)}}{2a({y}^{T},x)},\end{eqnarray} \tag{ 13 }$

where ${\rm{\Delta }}{({y}^{T},x)=b({y}^{T},x)}^{2}-4a({y}^{T},x)c({y}^{T},x)$ . Note that for each pair of right and left eigenvectors we have two possible solutions, but just one of them is an eigenvalue of ${\bf{Q}}(\lambda )$ . Furthermore, differentiating equation (12) by p we obtain information on how the eigenvalues change as a function of p. Formally we have

$\begin{eqnarray}&&\displaystyle \frac{\partial {y}^{T}{\bf{Q}}(\lambda )x}{\partial p}={y}^{T}\displaystyle \frac{\partial {\bf{Q}}(\lambda )}{\partial p}x+{y}^{T}{\bf{Q}}(\lambda )\displaystyle \frac{{\rm{d}}{x}}{{\rm{d}}{p}}+\displaystyle \frac{{\rm{d}}{y}}{{\rm{d}}{p}}{\bf{Q}}(\lambda )x=0,\end{eqnarray} \tag{ 14 }$

where

$\begin{eqnarray}&&\displaystyle \frac{\partial {\bf{Q}}(\lambda )}{\partial p}=2\lambda {\bf{A}}\displaystyle \frac{{\rm{d}}\lambda }{{\rm{d}}{p}}+\displaystyle \frac{{\rm{d}}\lambda }{{\rm{d}}{p}}{\bf{B}}+\lambda \displaystyle \frac{\partial {\bf{B}}}{\partial p}+\displaystyle \frac{\partial {\bf{C}}}{\partial p}.\end{eqnarray} \tag{ 15 }$

Note that the eigenvalues and eigenvectors are also a function of p. For continuity, two different eigenvalues may cross each other when varying p. Observe that for non-crossing points the relations $\displaystyle \frac{{{\rm{d}}{y}}^{T}}{{\rm{d}}{p}}{\bf{Q}}(\lambda )x\,=\,0$ and ${y}^{T}{\bf{Q}}(\lambda )\displaystyle \frac{{\rm{d}}{x}}{{\rm{d}}{p}}=0$ holds, since the derivatives are bounded for non-crossing points. Note that on the eigenvalue crossings we have two eigenvectors associated to the same eigenvalue. Thus, two solutions for the derivatives. Next, isolating the derivative of λ we have

$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\lambda }{{\rm{d}}{p}}=\displaystyle \frac{{y}^{T}\left(-\lambda \tfrac{\partial {\bf{B}}}{\partial p}-\tfrac{\partial {\bf{C}}}{\partial p}\right)x}{{y}^{T}(2\lambda {\bf{A}}+{\bf{B}})x}.\end{eqnarray} \tag{ 16 }$

In practical terms, this relation can be applied to drive a system through different structural or dynamical regimes. For example, considering the adjacency matrix, one can use this equation in order to chose an edge or set of edges to be removed (or weighted) in order to optimally reduce or increase the leading eigenvalue. As a consequence, the critical point of spreading processes, such as epidemic spreading, would also change. Obviously the matrix under study depends on the process. Another application is to design a numerical method to follow the correct eigenvalues as a function of p in a problem that might present eigenvalues crossings [7, 18].

3.2. Bounds

If x is an eigenvector (left or right), then ${x}^{T}{\bf{Q}}(\lambda )x\,=\,0$ holds and we can use it to find bounds to equation (13). Multiplying ${\bf{Q}}(\lambda )$ by ${\bf{D}}$ , we get a monic polynomial, then if we bound $b({x}^{T},x)$ and ${\rm{\Delta }}({x}^{T},x)$ , both solutions of equation (13) will also be bounded. Thus, to bound those terms, we can use the numerical range of the matrices to which they are related. The numerical range is formally defined for any matrix ${\bf{X}}$ as $F({\bf{X}})=\{{x}^{T}{\bf{X}}x:x\in {\mathbb{C}}\mathrm{and}\ {x}^{T}x=1\}$ . Additionally, $\sigma ({\bf{X}})\subseteq F({\bf{X}})$ , where $\sigma ({\bf{X}})$ is in the set of eigenvalues of ${\bf{X}}$ . Additionally, if ${\bf{X}}$ is an Hermitian matrix ${x}^{T}{\bf{X}}x$ is the Rayleigh quotient of ${\bf{X}}$ , which implies ${\lambda }_{1}({\bf{X}})\leqslant {x}^{T}{\bf{X}}x\leqslant {\lambda }_{N}({\bf{X}})$ . Finally, to bound non-Hermitian matrices we use the relation of the spectral norm and the numerical range, given as $\tfrac{1}{2}{| | | {\bf{X}}| | | }_{2}\leqslant r({\bf{X}})\leqslant {| | | {\bf{X}}| | | }_{2}$ , where $r({\bf{X}})$ is its numerical radius, defined as $r({\bf{X}})=\mathop{\max }\limits_{{\parallel x\parallel }_{2}\,=\,1}| {x}^{* }{\bf{X}}x| \,=\max \{| z| :z\in F({\bf{X}})\}$ .

First, the term b(x^T, x) is bounded by

$\begin{eqnarray}&&-{| | | {\bf{B}}| | | }_{2}\leqslant b({x}^{T},x)\leqslant {| | | {\bf{B}}| | | }_{2},\end{eqnarray} \tag{ 17 }$

however, in many cases, ${\bf{B}}$ is an Hermitian matrix, which allow us to improve the bound as

$\begin{eqnarray}&&{\lambda }_{\min }({\bf{B}})\leqslant b({x}^{T},x)\leqslant {\lambda }_{\max }({\bf{B}})\end{eqnarray} \tag{ 18 }$

Interestingly, we remark that ${\bf{B}}$ is often related to the aggregated network, connecting both scales of the same structure.

Next, we evaluate ${\rm{\Delta }}{({x}^{T},x)=({x}^{T}{\bf{B}}x)}^{2}-4{x}^{T}{\bf{C}}x$ . Firstly, we analyze the term ${({x}^{T}{\bf{B}}x)}^{2}$ , by observing that: (a) $\min \{{\mu }_{i}\}\leqslant {x}^{T}{\bf{B}}x\leqslant \max \{{\mu }_{i}\}$ , (b) $\min \{{\mu }_{i}^{2}\}\leqslant x{{\bf{B}}}^{2}x\leqslant \max \{{\mu }_{i}^{2}\}$ and (c) $\min {\{| {\mu }_{i}| \}}^{2}\leqslant {({x}^{T}{\bf{B}}x)}^{2}\,\leqslant \max {\{| {\mu }_{i}| \}}^{2}$ , since $\min \{{\mu }_{i}^{2}\}=\min {\{| {\mu }_{i}| \}}^{2}$ , hence, from (b) and (c), bounding ${({x}^{T}{\bf{B}}x)}^{2}$ is equivalent to bound $x{{\bf{B}}}^{2}x$ . Secondly, we can factorize ${\rm{\Delta }}({x}^{T},x)={x}^{T}({{\bf{B}}}^{2}-4{\bf{C}})x$ and defining the matrix ${\rm{\Delta }}={{\bf{B}}}^{2}-4{\bf{C}}$ , we can focus on the problem ${x}^{T}{\rm{\Delta }}x$ instead of the initial definition of ${\rm{\Delta }}({x}^{T},x)$ , since both have the same bounds. In addition, since frequently we are interested in undirected networks and consequently symmetric matrices, we can also constraint ${\rm{\Delta }}({x}^{T},x)\geqslant 0$ since we already know that the spectra is real in this case. Therefore, we have

$\begin{eqnarray}&&0\leqslant {\rm{\Delta }}({x}^{T},x)\leqslant {| | | {\rm{\Delta }}| | | }_{2}.\end{eqnarray} \tag{ 19 }$

Observe that those bounds can be further improved when applied to the analysis of particular matrices (supra-adjacency, supra-Laplacian, and probability transition) since their particularities also impose constraints on the solutions and could be explored to improve the bounds.

3.3. Comments on symmetric problems: HQEP

As previously mentioned, if ${\bf{A}},{\bf{B}}$ and ${\bf{C}}$ are Hermitian, we have a special class of problems called HQEP [16]. The HQEP has interesting properties, for instance, if x is a right eigenvector associated with the eigenvalue λ, then it is also a left eigenvector of the same eigenvalue [16]. In order to take advantage of those properties one can interpret the original problem as an HQEP plus asymmetric perturbation. Thus, the matrix polynomial defined by the matrix coefficients in equation (9) is, in general, asymmetric. However, a class of problems that arise naturally is defined by ${{\bf{M}}}_{12}=p{\bf{I}}$ . In this case, the matrices ${\bf{A}}$ and ${\bf{B}}$ are Hermitian. Observe that ${\bf{C}}$ still might be asymmetric. Nonetheless, we can use the Toeplitz decomposition [19] in order to analyze a symmetric and simplified problem. This decomposition states that any square matrix can be uniquely written as the sum of an Hermitian ( ${\bf{X}}={{\bf{X}}}^{* }$ ) and a skew Hermitian matrix ( ${\bf{X}}=-{{\bf{X}}}^{* }$ ) as ${\bf{X}}=\tfrac{1}{2}({{\bf{X}}}_{{\bf{1}}}+{{\bf{X}}}_{{\bf{2}}}^{* })+\tfrac{1}{2}({{\bf{X}}}_{{\bf{1}}}-{{\bf{X}}}_{{\bf{2}}}^{* })$ . This allows us to decompose $p{\bf{C}}=\tfrac{1}{2}({{\bf{M}}}_{11}{{\bf{M}}}_{22}+{{\bf{M}}}_{22}{{\bf{M}}}_{11})+\tfrac{1}{2}({{\bf{M}}}_{11}{{\bf{M}}}_{22}-{{\bf{M}}}_{22}{{\bf{M}}}_{11})+{p}^{2}{\bf{I}}$ . In this way, we can re-write our original QEP into two parts, a HQEP, composed only by Hermitian matrices and a skew Hermitian matrix, that can be interpreted as a perturbation. As a consequence from the perturbation theory, the matrix $p{\bf{C}}$ of the HQEP is perturbed by $\tfrac{1}{2}({{\bf{M}}}_{11}{{\bf{M}}}_{22}-{{\bf{M}}}_{22}{{\bf{M}}}_{11})$ and such matrix norm goes to zero as the layers are more similar. From the Bauer and Fike theorem [19] we can write a quality function for the approximation of the perturbed matrix ${\bf{C}}$ as

$\begin{eqnarray}&&| \lambda -\hat{\lambda }| \leqslant \kappa ({\bf{U}})\left|\left|\left|\displaystyle \frac{1}{2}({{\bf{M}}}_{11}{{\bf{M}}}_{22}-{{\bf{M}}}_{22}{{\bf{M}}}_{11})\right|\right|\right|,\end{eqnarray} \tag{ 20 }$

where $\hat{\lambda }$ is the eigenvalue of ${\bf{C}}={{\bf{C}}}_{{\bf{H}}}+{{\bf{C}}}_{{\bf{S}}}$ (Hermitian, ${{\bf{C}}}_{{\bf{H}}}$ , and skew-symmetric, ${{\bf{C}}}_{{\bf{S}}}$ ), ${{\bf{C}}}_{{\bf{H}}}={\bf{U}}{\rm{\Lambda }}{{\bf{U}}}^{-1}$ and $\kappa (\cdot )$ is the condition number with respect to the matrix norm $| | | \cdot | | |$ . Considering the spectral norm ${| | | \cdot | | | }_{2}$ we have $\kappa ({\bf{X}})=\left|\tfrac{{\sigma }_{\max }({\bf{X}})}{{\sigma }_{\min }({\bf{X}})}\right|$ . If $\kappa ({\bf{U}})$ is near 1, small perturbations imply small changes on the eigenvalues. Conversely, large values of $\kappa ({\bf{U}})$ suggest a poor approximation. Note that, although this analysis concerns only the matrix ${\bf{C}}$ and not the whole QEP, it can be an estimate of the quality of the approximation.

In addition to the HQEP properties, the perturbation analysis also emphasizes an important multiplex property. We must note that the more similar the layers are, the closer to zero the norm $\left|\left|\left|\tfrac{1}{2}({{\bf{M}}}_{11}{{\bf{M}}}_{22}-{{\bf{M}}}_{22}{{\bf{M}}}_{11})\right|\right|\right|$ is. Moreover, we also have another criteria which is based on the commutativity of the matrices ${{\bf{M}}}_{11}$ and ${{\bf{M}}}_{22}$ . Moreover, observe the role of correlations in this approximation. If both layers are identical, they are obviously correlated and the problem is symmetric.

3.4. Singular $D={M}_{12}={M}_{21}$ and the limits for sparse inter-layer coupling

So far we have assumed a node-aligned multiplex, i.e., a multiplex network in which each node has a counterpart on every layer[2], fulfilling the invertibility of ${{\bf{M}}}_{12}$ , which is a necessary condition to formally define the problem. Nevertheless, we can use the limit of ${{\bf{D}}}_{{ii}}=\epsilon \to 0$ to obtain an approximation of the sparse coupling. Thus, note that equation (5) can be analyzed in two different steps, firstly calculating the limit of decoupled edges and then the rest of the system. The first limit is analyzed as follows. From 5 the absent edges are factorized as

$\begin{eqnarray}&&\begin{array}{l}{p}^{-1}{\epsilon }^{-1}\tilde{{\bf{D}}}(i){\lambda }^{2}-{p}^{-1}{\epsilon }^{-1}({{\bf{M}}}_{11}\tilde{{\bf{D}}}(i)+\tilde{{\bf{D}}}(i){{\bf{M}}}_{22})\lambda \\ \quad +\,{p}^{-1}{\epsilon }^{-1}{{\bf{M}}}_{11}\tilde{{\bf{D}}}(i){{\bf{M}}}_{22}-p{\bf{D}}=0,\end{array}\end{eqnarray} \tag{ 21 }$

where $\tilde{{\bf{D}}}(i)=\epsilon {{\bf{D}}}^{-1}$ . Multiplying equation (21) by $p\epsilon$ and using the following limit

$\begin{eqnarray}\mathop{\mathrm{lim}}\limits_{\epsilon \to 0}{[\epsilon {{\bf{D}}}^{-1}]}_{{jj}}=\left\{\begin{array}{ll}1 & \mathrm{if}\hspace{0.5cm}{[{{\bf{D}}}^{-1}]}_{{jj}}\in O(\epsilon ),\\ 0 & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 22 }$

we have

$\begin{eqnarray}&&\tilde{{\bf{D}}}{\lambda }^{2}-({{\bf{M}}}_{11}\tilde{{\bf{D}}}+\tilde{{\bf{D}}}{{\bf{M}}}_{22})\lambda +{{\bf{M}}}_{11}\tilde{{\bf{D}}}{{\bf{M}}}_{22}=0,\end{eqnarray} \tag{ 23 }$

where the term of order $p\epsilon {\bf{D}}$ vanishes in the limit of $\epsilon \to 0$ . Note that $\tilde{{\bf{D}}}={\mathrm{lim}}_{\epsilon \to 0}[\epsilon {{\bf{D}}}^{-1}]={\bf{I}}$ if both layers are decoupled. Therefore, the factorization of the polynomial equation yields to $({{\bf{M}}}_{11}-\lambda {\bf{I}})({{\bf{M}}}_{22}-\lambda {\bf{I}})=0$ , whose solutions are the union of the solution of the standard eigenvalue problem of each layer. An important consequence is that the number of eigenvalues that do not change as a function of p are twice the number of nodes that do not have a counterpart on other layer.

Equation (23) give us the solutions for nodes without a counterpart on the other layer. To calculate the remaining solutions we have to redefine the original problem in terms of the Moore–Penrose pseudoinverse, which is denoted by ${{\bf{X}}}^{\dagger }$ , for a general matrix ${\bf{X}}$ . Thus, denoting by $\bar{{\bf{D}}}={p}^{-1}{{\bf{D}}}^{\dagger }$ , we have ${\bar{{\bf{D}}}}_{{jj}}={p}^{-1}{{\bf{D}}}_{{jj}}^{-1}$ if ${{\bf{D}}}_{{jj}}\ne 0$ and ${\bar{{\bf{D}}}}_{{jj}}=0$ otherwise. Observe that the zeros of ${\bar{{\bf{D}}}}_{{jj}}$ are ones in ${\tilde{{\bf{D}}}}_{{jj}}$ . For simplicity, in the rest of the paper we assume that ${{\bf{M}}}_{12}$ is invertible, however, the strategy above presented can be applied in case it is not. From the computational viewpoint, the cost of calculating the whole spectra for a closed range of p might be reduced since we can separate it into two parts: (i) a constant subset and (ii) the remaining subset, which varies as a function of p.

4. Applications

We next apply our proposed formalism to study the supra-Laplacian and the supra-adjacency matrices. For the sake of completeness, let us explicitly define the supra-adjacency matrix in terms of its block matrices (adjacency matrix of the individual layers). Formally, the supra-adjacency matrix is defined as

$\begin{eqnarray}{\bf{A}}=\left[\begin{array}{cc}{{\bf{A}}}_{{\bf{a}}} & p{\bf{I}}\\ p{\bf{I}} & {{\bf{A}}}_{{\bf{b}}}\end{array}\right],\end{eqnarray} \tag{ 24 }$

where we weight differently the intra and inter-layer edges. The definition of the supra-Laplacian matrix is

$\begin{eqnarray}{\bf{L}}=\hat{{\bf{D}}}-{\bf{A}}=\left[\begin{array}{cc}p{\bf{I}}+{{\bf{L}}}_{{\bf{a}}} & -p{\bf{I}}\\ -p{\bf{I}} & p{\bf{I}}+{{\bf{L}}}_{{\bf{b}}}\end{array}\right],\end{eqnarray} \tag{ 25 }$

where $\hat{{\bf{D}}}$ is a diagonal matrix whose elements are ${{\bf{D}}}_{{ii}}=\sum {{\bf{A}}}_{{ij}}$ and the Laplacian matrices of the individual layers are denoted as ${{\bf{L}}}_{{\bf{a}}}$ and ${{\bf{L}}}_{{\bf{b}}}$ . Many dynamical processes are described using those matrices. For instance, synchronization of coupled oscillators and also diffusion processes, while the supra-adjacency matrix is intimately related to epidemic and information spreading. It is also noteworthy that many structural metrics are also directly extracted from the spectral properties of those matrices. For instance, the communicability, which can be easily written as a matrix function, or more specifically, as the exponential of the adjacency or supra-adjacency matrix. Here we are going to focus on the spectral properties of these matrices, and their behavior as a function of the coupling parameter p under different conditions.

In addition to the supra-adjacency and supra-Laplacian matrix, we also analyze the probability transition matrices in section 4.3, which can be used to describe classical random walks on networks. The analysis of such a matrix is left to the last section since it is mainly numerical. Note that the probability transition matrix has a well bounded spectra where $1={\lambda }_{1}\geqslant {\lambda }_{2}\geqslant {\lambda }_{3}\geqslant ...\geqslant {\lambda }_{N}\geqslant -1$ [20]. This characteristic imposes an extra challenge on the derivation of the bounds. Although we could not improve those bounds, we report an interesting spectral behavior found numerically.

4.1. Supra-Laplacian matrix

The simplest supra-Laplacian matrix can be built considering a diagonal coupling matrix ${M}_{12}=-p{\bf{I}}$ , where each node has a counterpart on the other layer and the coupling is homogeneous. This implies that the QEP is defined with the following coefficient matrices

$\begin{eqnarray}&&{\bf{A}}={\bf{I}},\end{eqnarray} \tag{ 26 }$

$\begin{eqnarray}&&{\bf{B}}=-({{\bf{L}}}_{{\bf{a}}}+{{\bf{L}}}_{{\bf{b}}}+2p{\bf{I}}),\end{eqnarray} \tag{ 27 }$

$\begin{eqnarray}&&{\bf{C}}={{\bf{L}}}_{{\bf{a}}}{{\bf{L}}}_{{\bf{b}}}+p({{\bf{L}}}_{{\bf{a}}}+{{\bf{L}}}_{{\bf{b}}}).\end{eqnarray} \tag{ 28 }$

It is noteworthy that the aggregated network, ${{\bf{L}}}^{+}={{\bf{L}}}_{{\bf{a}}}+{{\bf{L}}}_{{\bf{b}}}$ , appears naturally under this formalism. This is interesting since it is physically understandable. On the other hand, the term ${{\bf{L}}}_{{\bf{a}}}{{\bf{L}}}_{{\bf{b}}}$ , in the definition of ${\bf{C}}$ , is of not so direct interpretation. This system presents a structural transition, which can be directly derived from our formalism. This derivation is presented in section 4.1.1. Additionally, we can also obtain bounds for the spectra using the ideas discussed in section 3.2. Those improved bounds are derived in section 4.1.2, where we use the particular properties of a Laplacian matrix to improve our previous results. In section 4.1.3, we evaluate the spectra of the supra-Laplacian matrix as a function of p and also compare our previous results with sparse and heterogeneous couplings. Specifically, on the heterogeneous case we consider a coupling matrix ${M}_{12}=-p{\bf{D}}$ , where ${\bf{D}}$ is a diagonal matrix. The QEP of such matrix is defined by ${\bf{A}}={{\bf{D}}}^{-1},{\bf{B}}=-({{\bf{L}}}_{{\bf{a}}}{{\bf{D}}}^{-1}+{{\bf{D}}}^{-1}{{\bf{L}}}_{{\bf{b}}}+2p{\bf{I}})$ and ${\bf{C}}={{\bf{L}}}_{{\bf{a}}}{{\bf{D}}}^{-1}{{\bf{L}}}_{{\bf{b}}}+p({{\bf{L}}}_{{\bf{a}}}+{{\bf{L}}}_{{\bf{b}}})$ . The analysis of such QEP is not trivial, since the matrices are asymmetric. Instead, we can explore it numerically and compare with the homogeneous case, ${M}_{12}=-p{\bf{I}}$ .

4.1.1. Structural transitions

First of all, we discuss the structural transition on the Laplacian matrix, initially presented in [21]. Here we calculate the exact transition points using the QEP formulation, which can be easily derived. We remark that those transition points were also calculated in [18], however, using two different methods: (i) eigenvalue sensitivity analysis and (ii) the Shurs complement approach. Both derivations presented in [18] are quite complicated, unlike our approach, where the solutions are obtained using simple arguments. Note, instead, that our approach presents a different expression if compared to the expression in [18]. Although their equivalence was not mathematically proven, it was numerically verified.

To begin with, it is well known that $\lambda =2p$ is an eigenvalue of the supra-Laplacian. Moreover, the crossing points are a consequence of the crossing between this eigenvalue and the bounded part of the supra-Laplacian spectra, yielding to the so-called structural transitions. Thus, from our defined QEP, we have that

$\begin{eqnarray}&&\det ({\bf{Q}}(2p))=\det ({{\bf{L}}}_{{\bf{a}}}+{{\bf{L}}}_{{\bf{b}}})\det ({{\bf{L}}}_{{\bf{b}}}{{\bf{L}}}_{{\bf{a}}}{({{\bf{L}}}_{{\bf{a}}}+{{\bf{L}}}_{{\bf{b}}})}^{\dagger }-p{\bf{I}}),\end{eqnarray} \tag{ 29 }$

which has two possible solutions: (i) $\det ({{\bf{L}}}_{{\bf{a}}}+{{\bf{L}}}_{{\bf{b}}})=0$ , which is always true, since the sum of the Laplacian matrices of the layers is the Laplacian of the aggregated network, hence, it also has a null determinant and (ii) the solution of $\det ({{\bf{L}}}_{{\bf{b}}}{{\bf{L}}}_{{\bf{a}}}{({{\bf{L}}}_{{\bf{a}}}+{{\bf{L}}}_{{\bf{b}}})}^{\dagger }-p{\bf{I}})$ , which are the crossing points or eigenvalues of multiplicity larger than one. Since this expression is also an eigenvalue problem, but in terms of p, we can express the crossing points as ${p}^{* }={\lambda }_{i}({{\bf{L}}}_{{\bf{b}}}{{\bf{L}}}_{{\bf{a}}}{({{\bf{L}}}_{{\bf{a}}}+{{\bf{L}}}_{{\bf{b}}})}^{\dagger })$ . Note that there are N possible values of p that solve equation (29), each one representing one crossing. The first solution is trivial, at p = 0, the second is the so-called structural transition [21]. We also remark its relevance in some dynamical processes, as presented in [22]. As said before, this expression is different from the previous one presented in the literature, however both give the same result.

4.1.2. Bounds

Here, exploiting the ideas presented in section 3.2, we improve the former bounds using specific Laplacian properties, such as its semi-positiveness. Thus, the QEP of the supra-Laplacian can be bounded using the individual bounds of ${\bf{B}}$ , which is a semi-positive definite Hermitian matrix, leading to

$\begin{eqnarray}&&2p\leqslant -b({x}^{T},x)\leqslant 2p+{\lambda }_{\max }({{\bf{L}}}_{{\bf{a}}}+{{\bf{L}}}_{{\bf{b}}}).\end{eqnarray} \tag{ 30 }$

Additionally, the discriminant function can also bounded by

$\begin{eqnarray}&&\begin{array}{l}\min \{{x}^{T}({({{\bf{L}}}_{{\bf{a}}}-{{\bf{L}}}_{{\bf{b}}})}^{2}-2p({{\bf{L}}}_{{\bf{a}}}+{{\bf{L}}}_{{\bf{b}}})+4{p}^{2}{\bf{I}})x\}\leqslant {\rm{\Delta }}({x}^{T},x)\\ \quad {\rm{\Delta }}({x}^{T},x)\leqslant \max \{x({({{\bf{L}}}_{{\bf{a}}}-{{\bf{L}}}_{{\bf{b}}})}^{2})+4{p}^{2}{\bf{I}}\},\end{array}\end{eqnarray} \tag{ 31 }$

where the upper bound can be defined as a function of the spectral properties of ${({{\bf{L}}}_{{\bf{a}}}-{{\bf{L}}}_{{\bf{b}}})}^{2}$ . On the other hand, the lower bound can be further improved by realizing that the matrix ${\rm{\Delta }}={({{\bf{L}}}_{{\bf{a}}}-{{\bf{L}}}_{{\bf{b}}})}^{2}-2p({{\bf{L}}}_{{\bf{a}}}+{{\bf{L}}}_{{\bf{b}}})+4{p}^{2}{\bf{I}}$ , defined on section 3.2, is semi-positive definite for undirected networks, ${\rm{\Delta }}\succcurlyeq 0$ . In this way, ${({{\bf{L}}}_{{\bf{a}}}-{{\bf{L}}}_{{\bf{b}}})}^{2}-2p({{\bf{L}}}_{{\bf{a}}}+{{\bf{L}}}_{{\bf{b}}})\,+4{p}^{2}{\bf{I}}\succcurlyeq 0$ , hence ${({{\bf{L}}}_{{\bf{a}}}-{{\bf{L}}}_{{\bf{b}}})}^{2}\,+\,4{p}^{2}{\bf{I}}\succcurlyeq 2p({{\bf{L}}}_{{\bf{a}}}+{{\bf{L}}}_{{\bf{b}}})$ , implying that ${\lambda }_{i}({({{\bf{L}}}_{{\bf{a}}}-{{\bf{L}}}_{{\bf{b}}})}^{2}\,+\,4{p}^{2}{\bf{I}})\,\geqslant {\lambda }_{i}(2p({{\bf{L}}}_{{\bf{a}}}+{{\bf{L}}}_{{\bf{b}}}))$ ⁶ .From these properties, we can establish the lower bound as 4p². Formally,

$\begin{eqnarray}&&4{p}^{2}\leqslant {\rm{\Delta }}({x}^{T},x)\leqslant {\lambda }_{\max }({({{\bf{L}}}_{{\bf{a}}}-{{\bf{L}}}_{{\bf{b}}})}^{2})+4{p}^{2}.\end{eqnarray} \tag{ 32 }$

The previously obtained bounds imply that, in the asymptotic analysis, we have ${\rm{\Delta }}({x}^{T},x)\in {\rm{\Theta }}({p}^{2})$ . Moreover, note that lower and upper bounds converge to each other as the layers become similar. On the extreme case of identical layers we have ${\rm{\Delta }}({x}^{T},x)=4{p}^{2}$ . Finally, combining the previously obtained bounds we have,

$\begin{eqnarray}&&0\leqslant {\lambda }^{-}({x}^{T},x)\leqslant \displaystyle \frac{1}{2}{\lambda }_{\max }({{\bf{L}}}_{{\bf{a}}}+{{\bf{L}}}_{{\bf{b}}})\end{eqnarray} \tag{ 33 }$

and

$\begin{eqnarray}\begin{array}{rcl}2p & \leqslant & {\lambda }^{+}({x}^{T},x)\,\\ & \leqslant & p+\displaystyle \frac{{\lambda }_{\max }({{\bf{L}}}_{{\bf{a}}}+{{\bf{L}}}_{{\bf{b}}})}{2}+\displaystyle \frac{\sqrt{{\lambda }_{\max }({({{\bf{L}}}_{{\bf{a}}}-{{\bf{L}}}_{{\bf{b}}})}^{2})+4{p}^{2}}}{2}.\end{array}\end{eqnarray} \tag{ 34 }$

Particularly, we can analyze these bounds in terms of their asymptotic behavior (approximation). Thus, for a sufficiently large value of p they can be approximated to

$\begin{eqnarray}&&0\leqslant {\lambda }^{-}({x}^{T},x)\leqslant \displaystyle \frac{{\lambda }_{\max }({{\bf{L}}}_{{\bf{a}}}+{{\bf{L}}}_{{\bf{b}}})}{2},\end{eqnarray} \tag{ 35 }$

$\begin{eqnarray}&&2p\leqslant {\lambda }^{+}({x}^{T},x)\leqslant 2p+\displaystyle \frac{{\lambda }_{\max }({{\bf{L}}}_{{\bf{a}}}+{{\bf{L}}}_{{\bf{b}}})}{2}.\end{eqnarray} \tag{ 36 }$

Therefore, from the asymptotic point of view we have ${\lambda }^{-}(x)\in {\rm{\Theta }}(1)$ and ${\lambda }^{+}(x)\in {\rm{\Theta }}(p)$ .

Figure 1 is an example of the eigenvalues as a function of the coupling parameter p for a multiplex network composed by two Erdös–Renyi layers. On the first layer we have an average degree $\langle k\rangle =12$ , while the second has $\langle k\rangle =16$ and both layers have $n={10}^{3}$ nodes.

Figure 1. Refer to the following caption and surrounding text. — **Figure 1.** The supra-Laplacian eigenvalues, $\lambda ({\bf{L}})$ , varying the coupling parameter p for a multiplex network composed by two Erdös–Renyi layers with n = 10³ nodes. The first layer has an average degree of $\langle k\rangle =12$ and the second of $\langle k\rangle =16$ . The upper bounds and lower bounds are dashed and continuous lines, respectively.
Download figure:
Standard image High-resolution image

**Figure 1.** The supra-Laplacian eigenvalues, $\lambda ({\bf{L}})$ , varying the coupling parameter p for a multiplex network composed by two Erdös–Renyi layers with n = 10³ nodes. The first layer has an average degree of $\langle k\rangle =12$ and the second of $\langle k\rangle =16$ . The upper bounds and lower bounds are dashed and continuous lines, respectively.
Download figure:
Standard image High-resolution image

4.1.3. Spectral properties as a function of the coupling parameter p

In this section, we focus on the eigenvalues behavior as a function of the coupling parameter, ${\lambda }_{i}(p)$ . First of all, we apply the concepts of section 3.1 regarding the derivative of Q(λ). Consider the simplest case, where ${\bf{D}}={\bf{I}}$ . In this case, we have a monic polynomial matrix, where ${\bf{B}}$ depends on the aggregated network, which is semi-positive definite. Moreover, ${\bf{C}}$ is a matrix that contains the product of the Laplacian of both layers, accounting the similarities between them. In this way, equation (16) can be expressed as

$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\lambda }{{\rm{d}}{p}}=\displaystyle \frac{(2\lambda {y}^{T}x-\hat{b}({y}^{T},x))}{(2(\lambda -p){y}^{T}x-\hat{b}({y}^{T},x))},\end{eqnarray} \tag{ 37 }$

where $\hat{b}({y}^{T},x)={y}^{T}({{\bf{L}}}_{{\bf{a}}}+{{\bf{L}}}_{{\bf{b}}})x$ and ${y}^{T}x=\cos (\theta )$ is the cosine of the angle between left and right eigenvectors of our QEP. Observe that part of the spectra has $\displaystyle \frac{{\rm{d}}\lambda }{{\rm{d}}{p}}\to 0$ , while the other part has $\displaystyle \frac{{\rm{d}}\lambda }{{\rm{d}}{p}}\to 2$ as p increases, which can be proved as follows. First of all, let us suppose that λ is a constant function of p, then $\displaystyle \frac{{\rm{d}}\lambda }{{\rm{d}}{p}}\to 0$ since the denominator grows as a function of p and the numerator is bounded, by supposition. Second of all, let us suppose that λ grows with p^r, where r < 1. Thus, $\displaystyle \frac{{\rm{d}}\lambda }{{\rm{d}}{p}}\to 0$ , using the same arguments, since the linear function of the denominator dominates it. On the other hand, if r = 1 we have $\displaystyle \frac{{\rm{d}}\lambda }{{\rm{d}}{p}}\to 2$ , since both, the numerator and the denominator, grow linearly. Finally, with p^r, where r > 1, both the denominator and numerator are dominated by p^r, implying that, for large p, the derivative $\displaystyle \frac{{\rm{d}}\lambda }{{\rm{d}}{p}}\to 1$ , which is also a contradiction, as it was supposed to be a linear function of p. Thus, it implies that the derivatives of λ, for large values of p, cannot grow faster than linearly and their growth will be one of two values, 0 or 2. Interestingly, note that these results are in alignment with the previously obtained bounds. In addition, such behavior is exemplified in figure 1.

Although for the simplified case we have two possible solutions at large p, observe that the above arguments fail for the case of general coupling matrix. From equation (16) and the definition of the Laplacian QEP we conclude that only the denominator of equation (16) changes for a different choice of ${\bf{D}}$ , since the terms that have dependencies on ${\bf{D}}$ vanish in the partial derivatives of the numerator. The denominator follows the general form ${y}^{T}(2\lambda {{\bf{D}}}^{-1}-2p{\bf{I}}-{{\bf{L}}}_{{\bf{a}}}{{\bf{D}}}^{-1}-{{\bf{D}}}^{-1}{{\bf{L}}}_{{\bf{b}}})x$ . In this way, different coupling weights can change the behavior of each eigenvalue differently for large p. For instance, if ${\bf{D}}={\bf{I}}$ the spectral distribution for large p is bimodal, however, if ${\bf{D}}=p\left(\tfrac{n}{{\sum }_{i}^{n}i}\mathrm{diag}(1,2,\,\ldots ,\,n)\right)$ ⁷ , this behavior changes completely and the eigenvalues change with different rates, presenting a 'continuous' bulk. This argument is valid for infinity size networks, since for finite size networks for a large p gaps between eigenvalues may appear due to different rates of growth (as a function of p). An example of this is shown in figure 2. Furthermore, we also found an empirical function that seems to bound the spectra as a function of p in this experiment. The lower bound is trivial, since it is a semi-positive definite matrix. The upper bound can be obtained correcting $\tilde{p}=\max \{\mathrm{diag}({\bf{D}})\}p$ , hence

$\begin{eqnarray*}&&0\leqslant \lambda \leqslant \tilde{p}+\displaystyle \frac{{\lambda }_{\max }({{\bf{L}}}_{{\bf{a}}}+{{\bf{L}}}_{{\bf{b}}})+\sqrt{{\lambda }_{\max }({({{\bf{L}}}_{{\bf{a}}}-{{\bf{L}}}_{{\bf{b}}})}^{2})+4{\tilde{p}}^{2}}}{2}.\end{eqnarray*}$

From figure 2 we observe that such bound is not as close to the largest eigenvalue as the homogeneous case.

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** The supra-Laplacian eigenvalues, $\lambda ({\bf{L}})$ , varying the coupling parameter p, considering the coupling matrix ${\bf{D}}=p\left(\tfrac{n}{{\sum }_{i}^{n}i}\mathrm{diag}(1,2,\,\ldots ,\,n)\right)$ —heterogeneous coupling. The dashed line is the adapted upper bound. The network considered in this figure has the same layers as used in figure 1.
Download figure:
Standard image High-resolution image

Figure 3. Refer to the following caption and surrounding text. — **Figure 3.** The supra-adjacency eigenvalues, $\lambda ({\bf{A}})$ , varying the coupling parameter p. The upper bounds and lower bounds are dashed and continuous lines, respectively. The network considered is the same as in figure 1.
Download figure:
Standard image High-resolution image

**Figure 3.** The supra-adjacency eigenvalues, $\lambda ({\bf{A}})$ , varying the coupling parameter p. The upper bounds and lower bounds are dashed and continuous lines, respectively. The network considered is the same as in figure 1.
Download figure:
Standard image High-resolution image

In addition to a non-homogeneous coupling matrix, the last case studied is the sparse coupling. The analytical part of this study was presented in section 3.4. As predicted, each uncoupled node implies a pair of eigenvalues that does not vary as function of p. Due to the nature of the Laplacian matrix, where just half of the eigenvalues vary with p, while the other half remain bounded, the set of bounded eigenvalues increases by one. For instance, if we have $\tilde{n}$ uncoupled nodes, the bounded part have $n+\tilde{n}$ eigenvalues, while the 'unbounded' part have $n-\tilde{n}$ eigenvalues. Note that the upper bound for the bounded set of eigenvalues is not $\tfrac{1}{2}{\lambda }_{\max }({{\bf{L}}}_{{\bf{a}}}+{{\bf{L}}}_{{\bf{b}}})$ anymore. Nevertheless, the general upper bound for ${\bf{D}}={\bf{I}}$ seems to be also an upper bound for the sparse problem, as numerically verified. The figures for these experiments were not shown as they are visually similar to figure 1.

4.2. Supra-adjacency matrix

Similarly to the supra-Laplacian case, here we also begin with the simplest case, i.e., the diagonal homogeneous coupling, and increase the level of complexity considering heterogeneous inter-layer weights and sparsity. Thus, in the simplest case we have ${M}_{12}=p{\bf{I}}$ , therefore, the QEP, equation (5), is defined by the following coefficient matrices

$\begin{eqnarray}&&{\bf{A}}={\bf{I}},\end{eqnarray} \tag{ 38 }$

$\begin{eqnarray}&&{\bf{B}}=-({{\bf{A}}}_{{\bf{a}}}+{{\bf{A}}}_{{\bf{b}}}),\end{eqnarray} \tag{ 39 }$

$\begin{eqnarray}&&{\bf{C}}={{\bf{A}}}_{{\bf{a}}}{{\bf{A}}}_{{\bf{b}}}-{p}^{2}{\bf{I}}.\end{eqnarray} \tag{ 40 }$

Note that, in a way similar to the Laplacian, ${\bf{B}}$ is also defined in terms of the aggregated network. On the other hand, the physical interpretation of ${\bf{C}}$ is still difficult due to the product ${{\bf{A}}}_{{\bf{a}}}{{\bf{A}}}_{{\bf{b}}}$ .

In section 4.2.1 we improve the bounds proposed in section 3.2. Then, in section 4.2.2, we evaluate the spectral properties of the supra-adjacency matrix as a function of p in three different contexts: (i) diagonal homogeneous coupling, (ii) diagonal heterogeneous coupling and (iii) sparse diagonal homogeneous coupling. Note that, in order to analyze the heterogeneous coupling we must consider the general QEP with the following coefficient matrices ${\bf{A}}={{\bf{D}}}^{-1},{\bf{B}}=-({{\bf{A}}}_{{\bf{a}}}{{\bf{D}}}^{-1}+{{\bf{D}}}^{-1}{{\bf{A}}}_{{\bf{b}}})$ and ${\bf{C}}={{\bf{A}}}_{{\bf{a}}}{{\bf{D}}}^{-1}{{\bf{A}}}_{{\bf{b}}}-{p}^{2}{\bf{D}}$ .

4.2.1. Bounds

Similarly to the analysis performed for the supra-Laplacian, here we also extend the ideas presented in section 3.2 to the supra-adjacency matrix. First of all, regarding the diagonal heterogeneous coupling case, ${\bf{D}}={\bf{I}}$ , we can also find the spectral distribution of the adjacency matrix. Beginning with matrix ${\bf{B}}$ , we can bound it using its eigenvalues as

$\begin{eqnarray}&&{\lambda }_{\min }({{\bf{A}}}_{{\bf{a}}}+{{\bf{A}}}_{{\bf{b}}})\leqslant -b(x)\leqslant {\lambda }_{\max }({{\bf{A}}}_{{\bf{a}}}+{{\bf{A}}}_{{\bf{b}}}).\end{eqnarray} \tag{ 41 }$

Interestingly, those are the eigenvalues of the aggregated network, which have a clear physical meaning. Similarly, for the discriminant we have

$\begin{eqnarray}&&\begin{array}{l}{\lambda }_{\min }({({{\bf{A}}}_{{\bf{a}}}-{{\bf{A}}}_{{\bf{b}}})}^{2})\leqslant {\rm{\Delta }}{({x}^{T},x)}^{2}\\ {\rm{\Delta }}({x}^{T},x)\leqslant {\lambda }_{\max }({({{\bf{A}}}_{{\bf{a}}}-{{\bf{A}}}_{{\bf{b}}})}^{2}).\end{array}\end{eqnarray} \tag{ 42 }$

Finally, combining those bounds we can bound both solutions by

$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}\displaystyle \frac{1}{2}({\lambda }_{min}({{\bf{A}}}_{{\bf{a}}}+{{\bf{A}}}_{{\bf{b}}})-\sqrt{{\lambda }_{max}({({{\bf{A}}}_{{\bf{a}}}-{{\bf{A}}}_{{\bf{b}}})}^{2})+4{p}^{2}})\leqslant {\lambda }^{-}\\ \leqslant \,\displaystyle \frac{1}{2}({\lambda }_{max}({{\bf{A}}}_{{\bf{a}}}+{{\bf{A}}}_{{\bf{b}}})-\sqrt{{\lambda }_{max}({({{\bf{A}}}_{{\bf{a}}}-{{\bf{A}}}_{{\bf{b}}})}^{2})+4{p}^{2}})\end{array}\end{array}\end{eqnarray} \tag{ 43 }$

and

$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{1}{2}({\lambda }_{\min }({{\bf{A}}}_{{\bf{a}}}+{{\bf{A}}}_{{\bf{b}}})+\sqrt{{\lambda }_{\max }({({{\bf{A}}}_{{\bf{a}}}-{{\bf{A}}}_{{\bf{b}}})}^{2})+4{p}^{2}})\leqslant {\lambda }^{+}\\ \quad \leqslant \,\displaystyle \frac{1}{2}({\lambda }_{\max }({{\bf{A}}}_{{\bf{a}}}+{{\bf{A}}}_{{\bf{b}}})+\sqrt{{\lambda }_{\max }({({{\bf{A}}}_{{\bf{a}}}-{{\bf{A}}}_{{\bf{b}}})}^{2})+4{p}^{2}}),\end{array}\end{eqnarray} \tag{ 44 }$

which asymptotically converge (as an approximation) to

$\begin{eqnarray}&&p\pm \displaystyle \frac{{\lambda }_{\min }({{\bf{A}}}_{{\bf{a}}}+{{\bf{A}}}_{{\bf{b}}})}{2}\leqslant {\lambda }^{\pm }(x)\leqslant p\pm \displaystyle \frac{{\lambda }_{\max }({{\bf{A}}}_{{\bf{a}}}+{{\bf{A}}}_{{\bf{b}}})}{2}.\end{eqnarray} \tag{ 45 }$

In other words, we have a bimodal spectral distribution for the adjacency matrix, where half of the eigenvalues grow linearly with p, while the other half decrease with the same rate.

Figure 3 is an example of the eigenvalues of the supra-adjacency matrix and the calculated bounds as a function of the coupling parameter p for a multiplex network composed by two Erdös–Renyi layers. The network used in this example is the same as that used in the supra-Laplacian case.

4.2.2. Spectral properties as a function of the coupling parameter p

In its general form, the first derivative is given as

$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\lambda }{{\rm{d}}{p}}=\displaystyle \frac{2{{py}}^{T}{\bf{D}}x}{(2\lambda {y}^{T}{{\bf{D}}}^{-1}x+b({y}^{T},x))},\end{eqnarray} \tag{ 46 }$

where x and y^T are, respectively, the right and left eigenvectors associated with the eigenvalue λ. Firstly, regarding ${\bf{D}}={\bf{I}}$ and using a similar approach as previously applied to the supra-Laplacian case. Thus, suppose that λ is a constant function of p or a function of p^r with r < 1, however, it would give us $\displaystyle \frac{{\rm{d}}\lambda }{{\rm{d}}{p}}\sim p$ , yielding to a contradiction. Then, suppose that it is a linear function of p, which implies $\displaystyle \frac{{\rm{d}}\lambda }{{\rm{d}}{p}}\to \pm 1$ , depending on the sign of the linear coefficient. Finally, suppose that it is a function of p^r with r > 1 we obtain that $\displaystyle \frac{{\rm{d}}\lambda }{{\rm{d}}{p}}\to 0$ , since the denominator grows faster than the numerator, which again is a contradiction. In this way, based on such analysis we conclude that the first derivative of λ can assume only $\displaystyle \frac{{\rm{d}}\lambda }{{\rm{d}}{p}}\to \pm 1$ .

Secondly, for the general case observe that both, the numerator and the denominator of equation (46) vary as a function of ${\bf{D}}$ . Additionally, ${\bf{D}}$ weights the product of the components of the eigenvectors, which allows the derivatives to assume more values, even a 'continuous bulk' instead of the bimodal distribution of the diagonal homogeneous case, similarly to the case discussed for the supra-Laplacian matrix. Here we also use the coupling matrix ${\bf{D}}=p\left(\tfrac{n}{{\sum }_{i}^{n}i}\mathrm{diag}(1,2,\,\ldots ,\,n)\right)$ . We show the spectral evolution as a function of p for the non-homogeneous case in figure 4. Similarly to the Laplacian case, there are evidences that the bounds can be corrected using $\tilde{p}=\max \{\mathrm{diag}({\bf{D}})\}p$ , hence

$\begin{eqnarray}&&\displaystyle \frac{1}{2}({\lambda }_{\min }({{\bf{A}}}_{{\bf{a}}}+{{\bf{A}}}_{{\bf{b}}})-\sqrt{{\lambda }_{\max }({({{\bf{A}}}_{{\bf{a}}}-{{\bf{A}}}_{{\bf{b}}})}^{2})+4{\tilde{p}}^{2}})\leqslant \lambda \end{eqnarray} \tag{ 47 }$

$\begin{eqnarray}&&\lambda \leqslant \displaystyle \frac{1}{2}({\lambda }_{\max }({{\bf{A}}}_{{\bf{a}}}+{{\bf{A}}}_{{\bf{b}}})+\sqrt{{\lambda }_{\max }({({{\bf{A}}}_{{\bf{a}}}-{{\bf{A}}}_{{\bf{b}}})}^{2})+4{\tilde{p}}^{2}}).\end{eqnarray} \tag{ 48 }$

Here we can also obtain a similar conclusion as for the supra-Laplacian case. From figure 4 we observe that the corrected bounds are not as close to the homogeneous case.

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** The supra-adjacency eigenvalues, $\lambda ({\bf{A}})$ , varying the coupling parameter p, considering the coupling matrix ${\bf{D}}=p\left(\tfrac{n}{{\sum }_{i}^{n}i}\mathrm{diag}(1,2,\,\ldots ,\,n)\right)$ , which is an heterogeneous coupling. The dashed and continuous red lines are the adapted upper and lower bounds, respectively. The network considered in this figure has the same layers as used in figure 1.
Download figure:
Standard image High-resolution image

Finally, we study the sparse coupling case, whose analytical evaluation was presented in section 3.4. As predicted, each uncoupled node imply in a pair of eigenvalues that does not vary as a function of p. In this way, for $\tilde{n}$ uncoupled nodes we have $2\tilde{n}$ eigenvalues in the central part of the spectra that do not vary as a function of p. Next, $n-\tilde{n}$ grows linearly with p, while the other $n-\tilde{n}$ eigenvalues with −p. This is shown in figure 5. Note that in this figure the horizontal lines bounding the central part of the spectra are not calculated, but numerically obtained and are only shown to serve as a reference.

Figure 5. Refer to the following caption and surrounding text. — **Figure 5.** The supra-adjacency eigenvalues, $\lambda ({\bf{A}})$ , varying the coupling parameter p, considering a sparse coupling matrix. The dashed and continuous red lines are the adapted upper and lower bounds, respectively. The dotted line was obtained numerically considering the largest value of p and are shown just as a reference. The network considered in this figure has the same layers as used in figure 1.
Download figure:
Standard image High-resolution image

4.3. Probability transition matrix

In this section, we evaluate the probability transition matrix, mainly focusing on its spectral properties as a function of the coupling parameter p. Due to the probabilistic nature of this matrix, we were not able to improve its bounds. Therefore, we mainly report numerical results.

Formally, the probability transition matrix is defined as

$\begin{eqnarray}\begin{array}{rcl}{\bf{P}}={\hat{{\bf{D}}}}^{-1}{\bf{A}} & = & \left[\begin{array}{cc}{({\tilde{{\bf{D}}}}^{A})}^{-1} & 0\\ 0 & {({\tilde{{\bf{D}}}}^{B})}^{-1}\end{array}\right]\left[\begin{array}{cc}{{\bf{A}}}_{{\bf{a}}} & p{\bf{I}}\\ p{\bf{I}} & {{\bf{A}}}_{{\bf{b}}}\end{array}\right]\\ & = & \left[\begin{array}{cc}{({\tilde{{\bf{D}}}}^{A})}^{-1}{{\bf{A}}}_{{\bf{a}}} & p{({\tilde{{\bf{D}}}}^{A})}^{-1}\\ p{({\tilde{{\bf{D}}}}^{B})}^{-1} & {({\tilde{{\bf{D}}}}^{B})}^{-1}{{\bf{A}}}_{{\bf{b}}}\end{array}\right],\end{array}\end{eqnarray} \tag{ 49 }$

where ${\tilde{{\bf{D}}}}_{{ii}}^{X}={k}_{i}^{X}+p$ , and X = { A, B} represents the label of each layer. It is known that this matrix models the classical random walk, where the walker chooses a neighbor based on the weights of its surrounding edges.

It is important to mention that in [23] the authors studied random walks on top of multiplex networks and analyzed them in terms of the normalized supra-Laplacian matrix. This matrix is defined as

$\begin{eqnarray}&&{{\bf{L}}}^{\mathrm{RW}}={\hat{{\bf{D}}}}^{-1}{\bf{L}}={\bf{I}}-{\hat{{\bf{D}}}}^{-1}{\bf{A}}={\bf{I}}-{\bf{P}},\end{eqnarray} \tag{ 50 }$

where ${\bf{L}}$ is the supra-Laplacian and ${\bf{P}}$ is the probability transition matrix. Note that the normalized supra-Laplacian matrix is intimately related to the probability transition matrix. In fact, their spectra are trivially related. Furthermore, we can also relate the spectra of the normalized Laplacian as follows

$\begin{eqnarray}&&{{\bf{L}}}^{\mathrm{Norm}}={\hat{{\bf{D}}}}^{-\tfrac{1}{2}}{\bf{L}}{\hat{{\bf{D}}}}^{-\tfrac{1}{2}}={\bf{I}}-{\hat{{\bf{D}}}}^{-\tfrac{1}{2}}{\bf{A}}{\hat{{\bf{D}}}}^{-\tfrac{1}{2}},\end{eqnarray} \tag{ 51 }$

where ${\bf{S}}={\hat{{\bf{D}}}}^{-\tfrac{1}{2}}{\bf{A}}{\hat{{\bf{D}}}}^{-\tfrac{1}{2}}$ has the same set of eigenvalues as ${\bf{P}}$ and if v is an eigenvector of ${\bf{S}}$ , then ${\hat{{\bf{D}}}}^{-1}v$ is an eigenvector of ${\bf{P}}$ associated with the same eigenvalue [20]. Note, however, that ${\bf{S}}$ is symmetric [20]. In the context of random walks in multiplex networks, in [13] the author used the normalized supra-Laplacian matrix. Here, in this section, we will study ${\bf{P}}$ , defined in equation (49).

Next, following our formalism, from equation (49), we can define our QEP in its monic form as

$\begin{eqnarray}&&{\bf{A}}={\bf{I}},\end{eqnarray} \tag{ 52 }$

$\begin{eqnarray}&&{\bf{B}}=-{({({\tilde{{\bf{D}}}}^{A})}^{-1}{{\bf{A}}}_{{\bf{a}}}+{{\bf{A}}}_{{\bf{b}}}{\tilde{{\bf{D}}}}^{B})}^{-1}),\end{eqnarray} \tag{ 53 }$

$\begin{eqnarray}&&{\bf{C}}={({\tilde{{\bf{D}}}}^{A})}^{-1}{{\bf{A}}}_{{\bf{a}}}{{\bf{A}}}_{{\bf{b}}}{({\tilde{{\bf{D}}}}^{B})}^{-1}-{({\tilde{{\bf{D}}}}^{A})}^{-1}{({\tilde{{\bf{D}}}}^{B})}^{-1}{p}^{2}.\end{eqnarray} \tag{ 54 }$

Note that such quadratic polynomial present some similarities with the one for the supra-adjacency matrix, however the probability transition matrix is not symmetric and the matrices ${\tilde{{\bf{D}}}}^{X}$ presents a dependency on p. This fact, associated with the natural bound for stochastic matrices, make the derivation of the spectral bounds more complicated than the previous cases. Here we focus on the spectral properties of the probability transition matrix as a function of the coupling parameter p.

4.3.1. Spectral properties as a function of the coupling parameter p

For the sake of completeness, let us study the spectral properties of the transition matrix as a function of the coupling strength p. This exercise is much more of an example than a practical application since we already know that the spectra are bounded on stochastic matrices, which does not allow unbounded grow. In figure 6 we present the spectra as a function of the coupling parameter p. The first observation is that, aside from being bounded, the growth rate of the eigenvalues is quite different from what was observed for the Laplacian and adjacency cases.

Figure 6. Refer to the following caption and surrounding text. — **Figure 6.** The probability transition matrix eigenvalues, $\lambda ({\bf{P}})$ , varying the coupling parameter p. The network considered is the same as that considered in figure 1.
Download figure:
Standard image High-resolution image

**Figure 6.** The probability transition matrix eigenvalues, $\lambda ({\bf{P}})$ , varying the coupling parameter p. The network considered is the same as that considered in figure 1.
Download figure:
Standard image High-resolution image

Firstly, lets proceed with the analysis of equation (16) aiming for an approximation, which qualitatively describes the $\lambda (p)$ . First of all, the partial derivative of ${\bf{B}}$ can be expressed as

$\begin{eqnarray}&&\displaystyle \frac{\partial {\bf{B}}}{\partial p}=-\displaystyle \frac{\partial }{\partial p}{({\tilde{{\bf{D}}}}^{A})}^{-1}{{\bf{A}}}_{{\bf{a}}}-\displaystyle \frac{\partial }{\partial p}{({\tilde{{\bf{D}}}}^{B})}^{-1}{{\bf{A}}}_{{\bf{b}}},\end{eqnarray} \tag{ 55 }$

where the term

$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial p}{({\tilde{{\bf{D}}}}^{X})}^{-1}=-{({k}_{i}^{X}+p)}^{-2}.\end{eqnarray} \tag{ 56 }$

Next, expanding the partial derivative of ${\bf{C}}$ we have

$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {\bf{C}}}{\partial p} & = & +\displaystyle \frac{\partial }{\partial p}{({\tilde{{\bf{D}}}}^{A})}^{-1}{{\bf{A}}}_{{\bf{a}}}{{\bf{A}}}_{{\bf{b}}}{({\tilde{{\bf{D}}}}^{B})}^{-1}\\ & & +{({\tilde{{\bf{D}}}}^{A})}^{-1}{{\bf{A}}}_{{\bf{a}}}{{\bf{A}}}_{{\bf{b}}}\displaystyle \frac{\partial }{\partial p}{({\tilde{{\bf{D}}}}^{B})}^{-1}\\ & & -\displaystyle \frac{\partial }{\partial p}{({\tilde{{\bf{D}}}}^{A})}^{-1}{p}^{2}{({\tilde{{\bf{D}}}}^{B})}^{-1}\\ & & -2p{({\tilde{{\bf{D}}}}^{A})}^{-1}{({\tilde{{\bf{D}}}}^{B})}^{-1}\\ & & -{({\tilde{{\bf{D}}}}^{A})}^{-1}{p}^{2}\displaystyle \frac{\partial }{\partial p}{({\tilde{{\bf{D}}}}^{B})}^{-1}.\end{array}\end{eqnarray} \tag{ 57 }$

All the expressions obtained so far are quite complicated to be analyzed in its exact form. Thus, we will proceed with an asymptotic analysis, aiming for a hypothesis of a possible formula that qualitatively describes the behavior of ${\lambda }_{i}$ as a function of p. In other words, we propose a formula that fits the expected asymptotic behavior, but we also expect it to work for smaller values of p. We must remark that this analysis is an approximation and, in order to verify its validity we perform numerical fittings and evaluate the obtained errors.

From the previously mentioned perspectives, the asymptotic behavior are ${\left(\tfrac{\partial {\bf{C}}}{\partial p}\right)}_{{ij}}\in O({p}^{-1}),{\left(\tfrac{\partial {\bf{B}}}{\partial p}\right)}_{{ij}}\in O({p}^{-2}),{({\bf{B}})}_{{ij}}\in O({p}^{-1})$ and, obviously, ${({\bf{A}})}_{{ij}}\in O(1)$ ⁸ . First of all, in figure 7 we present some examples of functions with different asymptotic behaviors. In (a) we present a function in O(1), showing that it can be bounded by a constant, in (b) we show some functions in O(p⁻¹), while in (c) two functions, one in O(p⁻²) and the other in O(p⁻³). Note that we can approximate some of the terms in equations (55)–(57) to the functions in figure 7. In figure 7, we also show examples of our guessed asymptotic behavior for

$\begin{eqnarray*}\begin{array}{rcl} & & {({\tilde{{\bf{D}}}}^{A})}^{-1}\sim \displaystyle \frac{1}{({k}_{x}+p)}\in O({p}^{-1})\\ & & 2p{({\tilde{{\bf{D}}}}^{A})}^{-1}{({\tilde{{\bf{D}}}}^{B})}^{-1}\sim \displaystyle \frac{p}{{({k}_{x}+p)}^{2}}\in O({p}^{-1})\\ & & {({\tilde{{\bf{D}}}}^{A})}^{-1}{p}^{2}\displaystyle \frac{\partial }{\partial p}{({\tilde{{\bf{D}}}}^{B})}^{-1}\sim \displaystyle \frac{{p}^{2}}{{({k}_{x}+p)}^{3}}\in O({p}^{-1}).\end{array}\end{eqnarray*}$

Next, in figure 7(c) we show examples for

$\begin{eqnarray*}\begin{array}{rcl} & & \displaystyle \frac{\partial }{\partial p}{({\tilde{{\bf{D}}}}^{X})}^{-1}\sim \displaystyle \frac{1}{{({k}_{x}+p)}^{2}}\in O({p}^{-2})\\ & & {({\tilde{{\bf{D}}}}^{A})}^{-1}{{\bf{A}}}_{{\bf{a}}}{{\bf{A}}}_{{\bf{b}}}\displaystyle \frac{\partial }{\partial p}{({\tilde{{\bf{D}}}}^{B})}^{-1}\sim \displaystyle \frac{1}{{({k}_{x}+p)}^{3}}\in O({p}^{-3}).\end{array}\end{eqnarray*}$

Note that we considered a single value of k_x, without considering products between different constants. Besides, just one term is considered. We remark that the main goal of this exercise is to have insights on the qualitative behavior of more complicated functions, such as equation (57). In other words, the performed approximations are not expected to quantitatively predict those terms, but qualitative represent and 'catch' the main behavior of those functions.

Figure 7. Refer to the following caption and surrounding text. — **Figure 7.** Example of functions with different asymptotic behaviors. Functions in O(1) in (a), O(p⁻¹) in (b) and in O(p⁻²) and O(p⁻³) in (c). In all plots we consider k_a=1 and k_b=100. Note that the asymptotic class does not change when we multiply the function by a constant.
Download figure:
Standard image High-resolution image

Firstly, recall that the spectra on stochastic matrices is bounded, which consequently restricts its derivatives. In other words, λ_i ∈ O(1). However, for the sake of the argument, let us suppose that ${\lambda }_{i}={c}_{1}{p}^{r}+O({p}^{r-1})$ , hence $\displaystyle \frac{{\rm{d}}{\lambda }_{i}}{{\rm{d}}{p}}={c}_{1}{{rp}}^{(r-1)}+O({p}^{r-2})$ , where r is an integer. Thus, comparing with equation (16), we have

$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\rm{d}}{\lambda }_{i}}{{\rm{d}}{p}} & = & {c}_{1}{{rp}}^{(r-1)}+O({p}^{r-2})\\ & = & \displaystyle \frac{(-2{c}_{1}{p}^{r}+O({p}^{r-1}))\times O({p}^{-2})+O({p}^{-1})}{{c}_{1}{p}^{r}+O({p}^{r-1})+O({p}^{-1})},\end{array}\end{eqnarray} \tag{ 58 }$

that can be rewritten as

$\begin{eqnarray}&&\begin{array}{l}({c}_{1}{{rp}}^{(r-1)}+O({p}^{r-2}))({c}_{1}{p}^{r}+O({p}^{r-1})+O({p}^{-1}))\\ \quad =\,(-{c}_{1}{p}^{r}+O({p}^{r-1}))\times O({p}^{-2})+O({p}^{-1}),\end{array}\end{eqnarray} \tag{ 59 }$

which simplifies to

$\begin{eqnarray}&&\begin{array}{l}{c}_{1}^{2}{{rp}}^{(2r-1)}+O({p}^{2r-2})+O({p}^{-1})\\ \quad =\,{c}_{2}{p}^{r-2}+O({p}^{r-3})+O({p}^{-1}),\end{array}\end{eqnarray} \tag{ 60 }$

which implies that $\displaystyle \frac{{\rm{d}}\lambda }{{\rm{d}}{p}}\to 0$ and r ≤ 0 since on the left-hand side we have a function in $O(\max \{2r-1,-1\})$ , while, on the right-hand side we have a function in $O(\max \{r-2,-1\})$ . This simple analysis suggests that $r\leqslant 0$ , for consistency. Note that we are not inferring anything regarding its 'velocity' (how fast it goes to zero). Such arguments reinforce that λ_i ∈ O(1), as previously mentioned. However, there are a huge class of functions that satisfies such restriction. In order to satisfy the so far established restrictions, let us suppose that

$\begin{eqnarray}&&\hat{{\lambda }_{i}}=\displaystyle \frac{{k}_{0}{p}^{2}}{{(p+{c}_{0})}^{2}}+\displaystyle \sum _{k=1}^{K}\displaystyle \frac{{\tilde{c}}_{k}{p}^{k-1}}{{(p+{c}_{k})}^{k}}=\displaystyle \frac{{k}_{0}{p}^{2}}{{(p+{c}_{0})}^{2}}+O({p}^{-1}),\end{eqnarray} \tag{ 61 }$

which is a function that satisfies our previous analysis. Thus, it also implies that

$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\hat{{\lambda }_{i}}}{{\rm{d}}{p}}=\displaystyle \frac{2{k}_{0}{c}_{0}p}{{(p+{c}_{0})}^{3}}-\displaystyle \sum _{k=1}^{K}\displaystyle \frac{{\tilde{c}}_{k}{p}^{k-2}({c}_{k}-{{kc}}_{k}+p)}{{(p+{c}_{k})}^{k+1}}\end{eqnarray} \tag{ 62 }$

$\begin{eqnarray}&&=\displaystyle \frac{2{k}_{0}{c}_{0}p}{{(p+{c}_{0})}^{3}}-\displaystyle \sum _{k=1}^{K}\displaystyle \frac{{\tilde{c}}_{k}{p}^{k-1}}{{(p+{c}_{k})}^{k+1}}+O({p}^{-3}),\end{eqnarray} \tag{ 63 }$

which yields $\displaystyle \frac{{\rm{d}}\hat{{\lambda }_{i}}}{{\rm{d}}{p}}\in O({p}^{-2})$ . Next, from equation (16) we have

$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\lambda }_{i}}{{\rm{d}}{p}}\sim \displaystyle \frac{O(1)\times O({p}^{-2})+O({p}^{-1})}{O(1)\times O(1)+O({p}^{-1})}=\displaystyle \frac{O({p}^{-1})}{O(1)}\in O({p}^{-1}).\end{eqnarray} \tag{ 64 }$

Note that it allows a set of possible solutions and, among them, it allows our initial supposition, equation (61), i.e., ${\lambda }_{i}\sim \hat{{\lambda }_{i}}$ . Note that $\hat{{\lambda }_{i}}\in O(1)$ and $\displaystyle \frac{{\rm{d}}\hat{{\lambda }_{i}}}{{\rm{d}}{p}}\in O({p}^{-2})$ , which is also in O(p⁻¹), as expected from equation (64). Besides, for the sake of visualization, in figure 7(a) we show two examples of the leading term of equation (61).

Next, we proceed with a numerical experiment, extracting some eigenvalues presented in figure 6 we perform a fitting aiming to obtain the same curve. We chose 5 eigenvalues: (i) the leading eigenvalue, λ₁ = 1, just as a reference and to emphasize that our proposed equation also works for that case, (ii) λ₃, the first eigenvalue on the bulk (note that there can be a crossing between λ₃ and λ₂, which would change the index of the eigenvalue—here we are not going to enter into details of this possible crossing behavior and, in order to avoid that, we chose to follow the third eigenvalue), (iii) λ_N, the smallest eigenvalue and (iv) the two intermediate eigenvalues ${\lambda }_{n-1}$ and ${\lambda }_{n+1}$ , where we just considered their values after the spectra divides into two parts. Note that the error of these two curves is expected to be larger than the previous cases. It is important to remark that, as previously mentioned, there can be a crossing between λ₃ and λ₂, but here we are looking at the main global behavior and such a change would not be a big source of error. In this way, we are showing that there is a set of parameters that approximates the spectra using the proposed equations. The proposed mentioned experiment does not serve as a proof, but it does serve as an evidence of such, or a similar, behavior.

Following the proposed pipeline, firstly, in order to obtain the fittings, we used the nonlinear least squares method, the Levenberg–Marquardt algorithm [24–26] and the least absolute residual robust regression. Additionally, all the initial conditions were set to one. In figure 8 we show the obtained fittings and the numerically obtained eigenvalues. Complementary, in table 1 we present the fitted parameters. Interestingly, we observe that the proposed approximation fits really well the observed curves, which can be objectively measured by means of the sum of squares due to error, whose values are also reported in table 1. Thus, the behavior of λ_i assumed in equation (61) seems to be a very good guess. Besides, we also observe that there seems to be a symmetry on the obtained parameters for λ₃ and λ_N, which are close. The only exception is k₀, since both have the same modulus, but with a different sign, as expected. The last important observation also regards the parameter $| {k}_{0}|$ . Note that such a parameter is very close to one on all the fittings, suggesting some underlying property of our formulation.

Figure 8. Refer to the following caption and surrounding text. — **Figure 8.** The probability transition matrix eigenvalues, λ_i(p), varying the coupling parameter p, for i = N and i = 3 (or i = 2 after the crossing, in order to have a continuous curve). The network considered is the same as that considered in figure 1. The dots are the obtained eigenvalues from eigendecomposition of ${\bf{P}}$ , while the continuous blue lines are fitted curves from equation (61), where we used just the first term on the summation, i.e., K = 1.
Download figure:
Standard image High-resolution image

**Figure 8.** The probability transition matrix eigenvalues, λ_i(p), varying the coupling parameter p, for i = N and i = 3 (or i = 2 after the crossing, in order to have a continuous curve). The network considered is the same as that considered in figure 1. The dots are the obtained eigenvalues from eigendecomposition of ${\bf{P}}$ , while the continuous blue lines are fitted curves from equation (61), where we used just the first term on the summation, i.e., K = 1.
Download figure:
Standard image High-resolution image

Table 1. Parameter values of the network reported in figure 8. The confidence intervals are given in parenthesis and the goodness of fit is measured by the sum of squares due to error (SSE).

Eigenvalue	k₀	c₀	c_k	${\tilde{c}}_{k}$	SSE
λ₁ = 1	1 (1, 1)	<10⁻⁵	<10⁻⁵	<10⁻⁵	<10⁻⁶
λ₃	0.9999 (0.9999, 1)	5.514 (5.502, 5.526)	8.618 (8.579, 8.657)	4.732 (4.712, 4.753)	< 10⁻⁴
λ_n−1	0.9985 (0.9975, 0.9996)	11.27 (11.18, 11.36)	−2.141 (−2.629, −1.653)	−0.1541 (−0.1839,−0.124 2)	< 10⁻²
λ_n+1	−0.9997 (−0.9999, −0.9995)	0.1273 (0.1186, 0.1361)	20.73 (20.61, 20.85)	23.58 (23.53, 23.63)	<10⁻⁴
λ_N	−1 (−1, −1)	5.344 (5.334, 5.353)	8.444 (8.412, 8.477)	−4.615 (−4.632, −4.598)	< 10⁻⁴

Finally, for the sake of completeness and for comparison reasons, we numerically evaluate the spectra of the probability transition matrix for the sparse and heterogeneous coupling cases. Regarding the sparsity, in figure 9, we present a similar experiment as done for the supra-Laplacian and supra-adjacency cases. Similarly to those experiments, here we also observe a group of eigenvalues that do not change as a function of p. Moreover, we also verified that for $\hat{n}=100$ decoupled nodes, we have $2\hat{n}=200$ eigenvalues that remain constant, validating the insights we obtained in section 3.4. Although the behavior observed is similar to the previously studied matrices, for the probability transition matrix we observe a slightly different behavior for intermediate values of p (here 1 < p < 10), where the intermediate eigenvalues change, forming the 'central bulk'.

Figure 9. Refer to the following caption and surrounding text. — **Figure 9.** The probability transition matrix eigenvalues, $\lambda ({\bf{P}})$ , varying the coupling parameter p. The coupling matrix is sparse. The network considered in this figure has the same layers as used in figure 1.
Download figure:
Standard image High-resolution image

Furthermore, we remark that for the heterogeneous coupling ( ${\bf{D}}=p\left(\tfrac{n}{{\sum }_{i}^{n}i}\mathrm{diag}(1,2,\,\ldots ,\,n)\right)$ ), if compared with the supra-Laplacian and supra-adjacency, a completely different behavior emerged. In the probability transition matrix case, the spectra seem to be always bimodal. This effect is shown in figure 10, where, for a large enough value of p, the eigenvalues tend to a constant. It is also noteworthy that the rate at which this phenomenon takes place is much slower than the rate of the homogeneous case, shown in figure 6.

Figure 10. Refer to the following caption and surrounding text. — **Figure 10.** The probability transition matrix eigenvalues, $\lambda ({\bf{P}})$ , varying the coupling parameter p, considering the coupling matrix ${\bf{D}}=p\left(\tfrac{n}{{\sum }_{i}^{n}i}\mathrm{diag}(1,2,\,\ldots ,\,n)\right)$ , i.e. an heterogeneous coupling. The network considered in this figure has the same layers as used in figure 1.
Download figure:
Standard image High-resolution image

5. Discussion and conclusions

From the developed theory, we applied and analyzed three different matrices: (i) the supra-Laplacian, (ii) the supra-adjacency and (iii) the probability transition matrix. In all these cases we have considered three different coupling schemes: (a) diagonal homogeneous coupling, ${\bf{D}}=p{\bf{I}}$ , (b) diagonal homogeneous sparse coupling and (c) diagonal heterogeneous coupling, ${\bf{D}}=p\left(\tfrac{n}{{\sum }_{i}^{n}i}\mathrm{diag}(1,2,\,\ldots ,\,n)\right)$ . Regarding the supra-Laplacian and the supra-adjacency matrices, on the first scenario, (a), we were able to extract some analytical results regarding the derivatives of the eigenvalues, which suggested a different behavior for the other two cases, (b) and (c). On the other hand, regarding the probability transition matrix, due to its stochastic nature, we were not able to go further with the analytical analysis. However, we followed an asymptotic analysis, proposing a function that describes the eigenvalues behavior. This function was validated with numerical fittings of the original spectra. Although it is just an approximation, it also helps us understand the nature of the phenomena behind this structure. Furthermore, we also reported the differences between the spectral distributions for large p, where we can have bimodal, multi-modal or even a continuous bulk for the adjacency and Laplacian cases, just changing the coupling matrices. On the sparse case, this analysis was analytically supported, while the other cases were explored numerically.

Our analysis pointed out some important features about multiplex systems. As a general observation, as we increase p we will find (roughly) three different structural phases, which might take place at different points for each structure and matrix. Thus, the structural phases of a multiplex network can be defined as: (i) decoupled phase, for small values of p, where the layers are virtually decoupled and act by themselves, with a neglectable interaction, (ii) multiplex/multilayer phase, where the system is coupled and the intra-layer edges play an important role and (iii) a network of layers phase, where the structure of the network of layers plays the major role. Note that, from the perturbation theory point of view, in the decoupled phase the eigenvalues are basically the union of the eigenvalues of the individual layers plus some perturbations. On the other extreme, in the network of layers phase, the intra-layer edges might be understood as the perturbation, since p ≫ 1. In this case, we can interpret the system as a set of n virtually disconnected small networks, whose structure is given by the network of layers (considering a multiplex case, where each node has a counterpart on the other layers). Finally, the most interesting scenario is the multiplex phase, where the inter and intra-layer topologies play a fundamental role on the dynamics.

Throughout our analysis, we were able to verify these regimes in the different matrices we evaluated. Although all of them showed this behavior, the differences between those matrices are also evident. Considering the supra-Laplacian and supra-adjacency matrices with diagonal homogeneous coupling, we observe that, for a large enough p, the spectral distribution is bimodal, while for the sparse case we have three bulk's, where the central one results from the nodes that do not have an inter-layer edge. Finally, on the diagonal heterogeneous coupling, we observe a completely different behavior, where the eigenvalues are distributed into a single bulk. We remark that, although we were not able to analytically quantify this last phenomenon, our analysis suggested such a behavior.

Furthermore, comparing those results with the ones obtained using the probability transition matrix, we observed a completely different behavior. On the diagonal, homogeneous or heterogeneous cases, the spectra seem to be bimodal for a sufficiently large p. Note that for the homogeneous case this convergence to the bulks is much faster than the heterogeneous case. This is an interesting phenomenon since it contrasts with the supra-adjacency and supra-Laplacian cases, where the heterogeneous coupling implies a 'continuous' bulk. It is noteworthy that our predictions for a central bulk for the uncoupled nodes are also fulfilled for the probability transition matrix.

The analysis performed here emphasize the importance of a proper study of the structural phases in different contexts. The sparse and heterogeneous cases might change completely the spectra (depending on p). Obviously, the analysis should also take into account the correct matrix since the structural changes are different from case to case. In other words, different matrices present their phases in different intervals (values of p). In this context, dynamical insights can also be useful for a better understanding.

Since all the analyzed matrices are also related to dynamical processes, the results reported here will directly impact on these processes too. Note that in [22] the authors analyzed diffusion processes in multiplex networks and found the so-called superdiffusion. This process is described by the supra-Laplacian matrix and it is intrinsically connected to the so-called structural transition of this matrix as pointed in [22] and latter discussed in [18], where the authors found the exact structural transition point. Furthermore, in [7], while studying epidemic spreading in multiplex networks, the authors verified this structural behavior in the analysis of the supra-adjacency matrix. Besides, it was also shown that it is intimately related to spreading processes and the layer-localization phenomena [7]. Thus, in the mentioned cases, the dynamical regimes can be understood as a consequence of the structural changes.

In summary, we have proposed a new mathematical formalism for the analysis of spectral properties in multiplex networks using the polynomial eigenvalue problem. This approach reduces the dimensionality of our matrices (coefficient matrices) at the cost of a higher order of the characteristic polynomial. This technique might seem counterintuitive at a first glance, but it reveals an underlying relationship between the eigenvalues of the matrices associated to multiplex structures. In contrast to single-layer networks, multiplex networks are defined as matrix functions since we are interested in weighting inter and intra-layer edges differently. Thus, they depend on the coupling parameter, here denoted by p. Therefore, it is of utmost importance to derive spectral bounds as a function of p and to obtain insights on their asymptotic behavior.

In addition to the so far discussed results, we also must comment the implications of our results to single-layer networks and also in other contexts. We remark that, aside from the multiplex networks, we can also have other types of regularities allowing the polynomial to be reduced and analytically evaluated. For instance, more general multilayer networks or some specific community structured single-layer networks. In fact, any problem where the structure can be divided into groups and whose groups have some regularities allowing the matrices ${\bf{B}}$ and ${\bf{C}}$ to be calculated. Note that one might even argue that the groups do not need to be at the same size since one can add phantom nodes in order to mathematically allow the groups to have the same size. Thus, from the results of section 3.4, the analysis would follow naturally. We remark that in this scenario the analysis should be carried out with an extra care on the physical meaning of the newly added eigenvalues due to the phantom nodes. Furthermore, in the community structured single-layer networks case note that some of our results are conditioned to the scenario where intra and inter community edges are weighted differently, playing the role of our coupling parameter p. We might also suppose that our formalism can be applied in dynamical contexts when we have two different dynamical sets of nodes. For instance, the chimera states in arrays of identical oscillators, where we have two groups of nodes: coherent and incoherent [27]. Finally, the numerically obtained results might also be of interest, giving insights in many different fields ranging from the structural analysis of single-layered community structured networks, to the study of dynamical processes.

We hope this work motivates the community to study in further details the structural behavior of multiplex and multilayer systems and their dynamical consequences. Besides, other matrices and processes might also be studied and evaluated. In this case, we believe that our formalism might also be helpful. Finally, we also hope to motivate studies on the analysis of eigenvectors, which is still lacking in the literature and are of great importance as they were shown to play a major role in dynamical processes [7, 8].

Acknowledgments

FAR acknowledge CNPq (grant 307974/2013-8) and Fapesp (grant 2013/26416-9) for financial support. GFA acknowledges Fapesp for the sponsorship provided (grants 2012/25219-2). YM acknowledges partial support from the Government of Aragón, Spain through grant E36-17R, and by MINECO and FEDER funds (grant FIS2017-87519-P).

: Appendix. Jordan triple

The definitions presented in this section are not studied in details here since our goal is focused on applications, however, we refer the reader to [14, 15] for more information on this class of problems. Here we reproduce some important definitions that might me useful for some readers, allowing them to extend the results presented in this work to other contexts.

Definition A.1. Jordan triple [15]: Denoting by ${\bf{J}}$ the Jordan matrix, $({\bf{X}},{\bf{J}},{\bf{Y}})$ is the Jordan triple of ${\bf{Q}}(\lambda )$ , where ${\bf{X}}$ and ${\bf{Y}}$ are the right and left eigenvectors. We also have that $({{\bf{Y}}}^{* },{{\bf{J}}}^{* },{{\bf{X}}}^{* })$ is the Jordan triple of ${\bf{Q}}{(\lambda )}^{* }$ . Note that the Jordan matrix, ${\bf{J}}$ , the diagonal blocks are the Jordan blocks. Besides, if all the eigenvalues are simple, ${\bf{J}}$ is a diagonal matrix, where ${{\bf{J}}}_{{ii}}={\lambda }_{i}$ . Additionally, observe that ${\bf{X}}\in {{\mathbb{R}}}^{n\times {nl}}$ , where ${\bf{X}}=[{x}_{1},\,\mathrm{...},{x}_{{nl}}]$ is composed by the right eigenvectors, while ${\bf{Y}}\in {{\mathbb{R}}}^{{nl}\times n}$ , where ${\bf{Y}}={[{y}_{1},...,{y}_{{nl}}]}^{T}$ is composed by its left eigenvectors. Finally, ${\bf{J}}\in {{\mathbb{R}}}^{{nl}\times {nl}}$ .

Definition A.2. The left set of eigenvectors ${\bf{Y}}$ [15]: The left eigenvectors can be defined in terms of the right eigenvectors and the matrix ${\bf{J}}$ as

$\begin{eqnarray}&&{\bf{Y}}={\left[\begin{array}{c}{\bf{X}}\\ {\bf{XJ}}\\ \vdots \\ {{\bf{XJ}}}^{l-1}\end{array}\right]}^{-1}\left[\begin{array}{c}0\\ \vdots \\ 0\\ {\bf{I}}\end{array}\right].\end{eqnarray} \tag{ A1 }$

Definition A.3. Useful relations for the order two polynomial eigenvalue problem: